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ABSTRACT 



The problem of tracking submarine targets with passive 
sonobuoys is mathematically modelled, A comprehensive study 
is made of all the information available in the acoustic 
signals picked ud by the sonobuoys and of the usefulness of 
this information in the estimation orocess. The presence of 
non 1 i nea r i t i es in the tracking model leads to the applica- 
tion of nonlinear estimation theory. Bayes formulation con- 
cepts are applied to generate approximate solutions and 
filtering algorithms^ and the well known extended Kalman 
filter equations and higher order filtering algorithms are 
obtained from this approach. 

The concept of partitioning the measurements is 
presented and shown to bring advantages in computing effi- 
ciency and also» for nonlinear measurement Sr in tracking 
accuracy. A graphical interoretation of the action of Kalman 
filters is developed and provides insight into the impor- 
tanceof each variable in the filtering process. 

Extensive simulations^ designed to test the performance 
of the algorithms developed# are presented in graphical form 
and analyzed. 
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SYMBOLS AND CONVENTIONS FOR TEXT 



1 - Underlined small letters represent vectors# thus# 

X denotes the vector x. 



2 - 


Capital letters are used to represent matrices. 


3 - 


A caret indicates an estimated value# e.g.# x 


denotes an 


estimateofx. 


a - 


Time points are represented by t, or# when 

k 



between parentheses# simply by k# e.g.# x(k) denotes the 



value of X 


at time t . 

k 


5 - 


Estimates specify the time of estimation and the 


amount of 


information used# e.g.# x(k!j) denotes the esti- 



mate of X at time t, taking into account all the measurements 

UD to and including t^ . When k = j only one time point is 

indicated# thus# x(k) aenotes the estimate of x at time t, 

~ — k 

taking into account all measurements up to and including tj^ . 



6 - 


Probability density functions are represented by 


the small 


letter p# i.e. denotes the conditional 


density of 


X given y. When no confusion is possible# this is 



simplified to p(x'y). 

7 - The expectation operator is represented by E [ ]; 

the variance by Varf 1. 



0 



£( k ) 

as(k) 
a b(k) 

a k ) 

b { k ) 

b^(k) 

b(k) 

value 

c 

d .(k) 

X 

e { k ) 
je ( k ) 

f (k) 

0 

f. (k) 

1 

G(k) 

H ( k ) 
m 

P(k) 



expected value of the estimated state vector at 



Sc 



random rate of chanoe of target speed. 



random rate of change of target heading. 



random rate of change of frequency emitted. 



heading of target at time tj^ . 



bearing measurement by buoy i. 



second moment of the state estimate around expected 
a ( k ) . 



average speed of sound, 
distance of target to buoy i. 

Kronecker delta; =lifk=j;=Oifk=j. 
estimation error vector. 

expected value of the estimation error. 

rest frequency emitted by the target. 

frequency measured by buoy i. 

gain matrix. 

observation matrix. 

compress i on/expans i on factor 

estimation error covariance matrix. 
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Q(k) 



state excitation covariance matrix 



R(k J 


measurement noise cavariance 


matrix. 


S ( k 1 


speed of target . 




s . (k) 

X 


relative velocity of target 


towards buoy i . 


“s'"’ 


standard deviation ofa (k). 

s 






standard deviation ofa (k). 

D 




Oj( k ) 


standard deviation ofa^ (k). 




T 


time delay. 




T . . 


time delay difference between the signals 


by buoys 


i and j . 




T 


time spent in prediction. 





rece i ved 



P 

T time soent in estimation, 

e 

V ( k ) vector of q random measurement noise signals. 

^(kJ vector of m random forcing inputs. 

jc ( k J state vector of dimension n. 

x(k) x-position of target. 

x.(k) x-position of buoy i. 

1 

£(k) estimated state vector. 

j^(kj In (n + 3 ) /2] -d i mens i ona 1 vector containing the state 
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variables and the distinct conponents of the matrix P 



y (k) 



y. (k) 
1 

£(k)‘ 

Zk 

i .e. 



y-position of target, 
y-position of buoy i. 



vector of q measurements. 

set of all measurements up to and including time t 
[z( 0 ), 2(1), ..., z(k)) . 
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INTRODUCTION 



This work was oriented to the problem of optimally 
estimating characteristics and/or parameters of a certain 
class of nonlinear# dynamic# stochastic systems from 
observed output sequences which are noise corrupted non- 
linear functions of the system states. 

In particular# attention was directed to the problem of 
real-time# short-range optimal localization and tracking of 
a submarine target by passive acoustical means. The most 
accurate and reliable estimates of the target's parameters 
(position# heading# speed) are sought by optimally process- 
ing the measurements obtained from the acoustic signals 
transmitted by the target itself and received by special 
sonobuoys . 

In the unclassified literature there are presently 
available methods of passive tracking of submarine targets. 
Reference tl) utilizes doool er-sh i f ted frequency measure- 
ments obtained from a group of sensors# (2) uses bearing and 
dopp 1 e r-sh i f ted frequency measurements obtained from one or 
more sonobuoys. These methods provide very good first 
approximations to the solution of the problem and many of 
their concepts and notation are used in this work. 



As a first step/ an attempt was made to produce a 



comprehensive study of all the information available in the 
acoustic signals picked up by the sonobuoys and its useful- 
ness in the estimation process. 

The estimation algorithms were developed taking into 
account the advantages of processing information as soon as 
it becomes available. Flexibility to move/ remove and 
include sensors and measurements at any time was obtained. 
Constraints were included to account for limitations of 
practical/ inexpensive and presently available sonobuoys and 
computing eguipment. 

In Chapter II the problem is described in mathematical 
terms using state space techniques and the initial assump- 
tions and constraints are given. The model used for the 
target is similar to the model presented in [21/ but with a 
different set of states. The measurements obtainea from the 
signals received by the sonobuoys are c h a r ac t e r i zed and 
their functional rel at i onsh i ps to the parameters of the tar- 
get are analyzed. 

Chapter III discusses the general theoretical solution 
for the problem and the difficulties of implementing it. 
Approximate solutions are then sought and processing equa- 
tions are developed. The special vectors/ matrices and 
relations c ha r ac t e r i s t i c of the tracking problem are 
prepared for the application of the filtering eouations. 



In Chapter IV the concept of partitioned measurements 



is introduced and analyzed. The advantages in accuracy^ com- 
puting efficiency and filter flexibility are stressed. 

Practical and graphical i nt erpret at i ons of the estima- 
tion process help in visualizing nonlinear problems and 
motivate the adoption of iterative techniques. Increased 
accuracy and possible convergence improvements are shown to 
result from this approach. 

Chapter V introduces the interactive computer program 
that allowed the application of the ideas described in the 
previous chapters to the tracking problem. Selected simula- 



t i ons are 


di scussed 


and presented 


in graphical 


f orm . 


This 


research 


covers the 


subject 


i n 


a much more 


comprehensive way 


than ear 1 i er 


works # 


m # 


12] # 13] . The 



results obtained by using the methods and ideas collected 
and developed in this study show many ways of improving the 
tracking accuracy^ of increasing computing efficiency# of 
reducing divergence problems and of obtaining faster conver- 
gence. 



The importance of having the buoys placed in adequate 
positions# of having an adequate model for the target# and 
of having frequent# varied and precise measurements is 
stressed. 
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The final chapter summarizes the results of this inves- 
tigation and presents the conclusions and suggestions for 
further study. 

- 

■ ■ 
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II. PROBLEM FORMULATION 



A, THE ESTIMATION PROBLEM 

A general nonlinear/ stochastic/ dynamic system observed 

✓ 

by a set of nonlinear/ non s t a t i ona r y measurements can be 
represented/ for our purposes/ in state-space discrete form 
by the equations 

x(k + l) = f(x(k)/w(k)/t, /t, ._) (2.1) 

— _ _ _ k k+1 

z.(k) = h (x_(k) / v(k) / tj^ ) (2.2) 

where x_( k ) is the state vector of dimension n/ 

w(k) is the vector of the m random forcing inputs/ 
^(k) is the measurement vector of dimension a/ 

^(k) is the vector of p random measurement noises/ 

J^( ) and h( ) are general vector functions. 

The general estimation problem consists of 

--- knowing the statistical c ha r ac t e r i s t i c s of w(k) and 

^(k) and having a statistical representation of the initial 

conditions ^(0) of the system/ 

--- processing the observations ^(0)/ z^( 1 ) / .../ _£ ( j ) 

in a statistically optimal way to obtain the best estimate/ 

in some sense/ of the system state x(k) at time t, . 

— k 

If k > J the processing is called prediction/ if k < J 
it is called smoothing and if k = j it is called filtering. 



In this work some departures from the completely general 
nonlinear case were taken. The first assumption is that the 
measurement noise is additive and that the system equations 
are linear with respect to the random input vector Equa- 
tions (2.1) and (2.2) can then be written as 



x(k + l) = f(x(k),t, ft.T ) + g(x(k)ft ft ).w(k) 

k Tc+1 ^ - k k+1 - 

(2.3) 

^(k + 1) = h(^(k)ftj^) + v(k) (2. a) 

where _f_( ) and h^( ) represent general vector functions and 

£( ) a general matrix of functions. 



Secondlyf the measurement noise and random forcing input 
are assumed to be uncorrelated in timef zero meanf discrete 
Gaussian sequences/ independent of one another and indepen- 
dent of the initial condition of the state of the system/ 

i .e . / 



E [w(k)J =0 
E (w(k)^''^( j )1 
E (w(k)^^( j )] 



E [v(k)J = 0 (2.5) 

Q(k). 6, . E [v(k)v'^( j )] = R(k), 6, 

kj - - kj 

0 Elx(0)w^(k)l = 0 E(v(k)x^(0)l = 0 



6 



kj 



i f k = j 




i f k ^ j 
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THE PASSIVE TRACKING PROBLEM 



Many practical situations can be formulated as estima- 
tion problems and characterized by the equations described 
above. The situation that motivated this study assumes a 
submarine target following (most of the time) an approxi- 
mately constant speed and constant heading path in a field 
of passive sonobuoys. 

The target unintentionally emits acoustic signals which 
are picked up by hydrophones. The transduced signals are 
sent up by wire to the buoys and then f requency-modu 1 ate VHP 
carriers which are transmitted by the buoys and received at 
a nearby ship or aircraft. After the recovery of the signals 
data processing equipment generates measurement values which 
contain numerical information about the target parameters. 

The first step in the analysis of the signals collected 
by the sonobuoys is the attempt to detect the presence of a 
real target. 

After a target is detected/ frequency spectrum informa- 
tion is added to intelligence data in an attempt to classify 
it. Other measurements and information from any other source 
are also used to obtain a first approximation to the 
target's parameters (position/ heading/ speed)/ generally 
through the use of least mean square techniques (11/ (3). 

The third phase of the process/ if a wartime condition 
doesn't exist/ is the tracking of the target. That iS/ the 
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determination of its path while reducing, if possible, the 
initial uncertainties about its parameters, in real-time, 
through the optimal processing of the sparse information 
provided by the sonobuoys. The passive characteristic causes 
severe limitations but is necessary to avoid revealing to 
the target the fact that it is being observed. 

This investigation was devoted to the tracking phase and 
it was assumed that the main characteristics of the buoys as 
well as their positions are known, 

C. THE MODEL OF THE TARGET 

The problem is, of course, three-dimensional. Neverthe- 
less the measurements are not a direct function of the 
target's deoth but of the difference in depth between the 
target and the hydrophones. Anticipating the accuracy of the 
measurements and the precision of the data processing eauip- 
ment and algorithms, it is justifiable to consider only two 
dimensions initially, adding depth as the third dimension 
whenever necessary ana with no conceptual aifficulty. 

The target's path is frequently one with nearly constant 
speed and heading, disturbed by currents and random thrust 
and control variations. Intentional maneuvers which have no 
evasive purposes are normally simple and smooth. 

A basic plant having as states two position coordinates 
(x and y), heading (b) and speed (s) of the target, describ- 
ing a constant velocity path, as shown in Fig 1, seems 
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adequat e 



Random forcing inputs should be considered to account 
for noise processes and imperfections of the model » such as 
target maneuvers. 

As suggested by (2J these random forcing functions can 

be approximately represented by two independent zero-mean/ 

pi ecewi se-constant random rates of changer a a, » act- 

s b 

ing on the speed and heading of the target. 



The formulation of this basic plant in discrete state- 
space formr similar to Equation (3.3)/ can be obtained in 
the following way 



rt 



s(k+l) - s(k) = 



k+1 



a ( t ) dt = ( t - t ) . a (k ) 
s k+1 k s 



ft 



b(k+l) - b(k) = 



k+1 



a (t) dt = (t, , - t ).a (k) 

b k+1 k b 



x(k+l) - x(k) = 

ft. 



k+1 



s(t) cos b(t) dt = 



'k+1 



(s(k) + ( t - t, )a (k)l cos(b(k) f 

k s 



+ (t - t,)a(k)l dt = 
k b 



'k+1 



s(k) cos b(k) dt + 



ft 



k+1 



•^k+1 



(t - t, )a (k) cos b(k) dt - 
k s 



(t - t, )a,(k) s(k) sin b(k) dt + 
k b 



bk 



■ ” T,) s(k) cos b(k) + 

k+1 k 
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Figure 1 : Geometry of target and sensors . 
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+ 4 ft, , - t. ( a (k) cos b(k) - 

2 k+1 k s 

- a (k) s(k) sin b(k)J + HOT 
b 

where HOT represents Higher-Order Terms. 



In the same way» 



y(k + l ) 



- y(k) 



'*^k+l 

s(t) sin b(t) dt = 



(t, , - t) s(k) sin b(k) + 

k+1 k 

+ (t, -t)^ [a(k)sinb(k) + 

2 k+1 k s 

+ “.(k) s(k) cos b(k)) + HOT 

D 



( 2 . 6 ) 



(2.7) 



K (t,,T - t, )a (k) and ( t, . , - t, ) a, (k) are suffi- 

k+1 k s k+1 k b 

ciently small so that the higher order terms in Equations 
(2.6) and (2.7) can be neglected^ one has» in vector form 



x(k + l) = f(x(k)^t ft, ) + g(x(k)ft, ft, ).w(k) 
— k k+1 k k+1 — 



where 



x(k) = 



” X ( k )^ 

y ( k ) 
s(k) 
b(k) 



( 2 . 8 ) 



w(k) = 



a (k) 
s 

ct^ (k) 
b 



(2.9) 
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( 2 . 10 ) 



f(x(k)ft 't, J ■ 
- - k+1 k 



x(k) ♦ (t - t ) s(k) cos b(k) 

k+1 k 

y(k) + (t - t ) s(k) sinb(k) 

k+1 k 

s(k) 

b(k) 



and 



g ( X ( k ) , t / 1 ) = 
— — k+1 k 



2 ^ ^ k+1 



— ( t 
2 ^ k+1 



^^k+1 " ^k^ 



( 2 . 11 ) 



b(k) 


- 




k+1 


- t, )^ s(k) 
k 


s 1 n 


b(k) 


b(k) 






k+1 


- t, )^ s ( k ) 
k 


cos 


b(k) 



k+1 ’ ^k^ 



When frequency measurements are to be processed 

directly^ as shown in Section IIfO,l»o» an additional state 

variable must be included --- the rest frequency f emitted 

o 

by the t a rget . 



This frequency is assumed to be approx i mat e 1 y constant 

with a small random disturbance that is assumed to be 

zero-mean^ p i ecew i se-const ant » and independent of a and a. . 

s b 

The expanded state variable formulation is given by 



x®(k) = 



X (k) 



f (k) 
o 



( 2 . 12 ) 
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w®(k) = 



w(k) 



a^( k ) 



(a 



X® (k + n = f® (x®(k) ,t, ) + 

— — — 



^ g ( k ) f t , , - / 1 ).w^(k) 

^ — k+1 k — 



where 



- f >L ( •< J ' t^+1 ' ^ ■ 



f(x(k)ft, , # t, ) 
k+1 k 



f (k) 
o 



g® (x®(k)/t, _ /t, ) 
— — k+1 k 



£(£( k ) » t 



k+1 



t, ) 



0 





0 



(t 



k+1 



D. MEASUREMENTS 

The songbuoys can perform one or both of the 
tasks: 

--- Dick up underwater sound signals 

indicate the approximate direction of the 
sum of the received acoustic signals. 



. 13 ) 



(2. la) 



( 2 . 15 ) 



( 2 . 16 ) 



foil ow i nq 



vectorial 
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The target emits a signal r(t) which is distorted and 
modulated by the propagation medium between the source and 
the hydrophone/ and by extraneous sound sources present in 
the ocean. 

When all of these disturbances are of relatively mild 
strength/ the signal received by a stationary buoy i will be 
approximately an attenuated and delayed version of the pre- 
viously emitted signal/ i.e./ 

r (t) = a (t).r(t - x (t)) (2.t7) 

i i i 

where a (t) is an attenuation factor and x. (t) is the time 
i 1 

delay gi ven by 

T (t) = d (t - X (t)) / c (2.18) 

i i i 

that iS/ the distance d between the target and the buoy at 

i 

the time of emission divided Dy the average velocity of 
sound in the medium. 

Since the target velocity is much smaller than the sound 
velocity/ for small distances and delays one has 

d.(t - X (t)) = d (t) + s (t).T (t) (2.19) 

i- i i i i 

where s (t) is the relative speed of the target toward buoy 
i 

i as shown in Fig 2 . 

Thus/ from Equations (2.18) and (2.19) one has 

x.(t) = d. (t) / (c - s (t)) (2.20) 

11 i 
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position of target 







Figure 2 : Target-sensors geometry - distances 
and relative velocities. 
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and the received signal is given by 



r.(t) = a.(t).r[t - d.(t)/(c - s,(t))J (2.21) 

1-1 1 1 

If this signal is recorded from t=t tot= t + 

o o 

one has recorded 





r . C t 
1 ■ o 


f 


t ’ ) 


= a ( t 
i o 


f t ' 


) . r It + t * - 

o 












•d ( t 
i 1 


+ 

0 


t')/(c - s (t + t 
i o 


• ) )] 




0 ^ t 


' < 


T 












If during 


this period 


T the relative 


target speed 


towards buoy 


i doesn't change 


significantly/ 






s ( t 
1 o 




t ' ) 


= s ( t ) 
i o 










d . ( t 
1 o 




t ’ ) 


= d (t ) 
i o 


- t 


' .s (t ) 
i o 




and 


hence 
















r . ( t 
1 • o 


+ 


t • ) 


= a ( t 
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i o 



In Equation (2.22)/ the second term inside the brackets 
is a fixed time delay. The third term is a variable time 
delay and can be seen as a comoress i on/exoans i on term, show- 
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ing the different value of the variable time for the origi- 
nal and the received signals. 

If the emitted signal had a single freguency component 
one would have 

r (t + t') = a (t + t').cos(w (t f t')) 
i • o i o i o 

where 

cos(w.(t t')) = cos(<(> w .t') = 

i o io i 

= cos fw .t - w ,d (t )/(c - s (t )) + 
o o o i o i o 

+ w .c.tV(c - s (t ))] = 
o i o 

= cos { (f) + w .c.t'/(c - s (t ))] 

io o i o 

and the received freguency would be 

w = w.c/(c-s) 
i o i 

f = f .c / (c - s ) (2.23) 

i o i 

The signals collected by various sonobuoys are sent to a 
ship or aircraft normally eouipped with equipment and skills 
to extract the information described below. 

I . I so 1 a t ed Buoys 

By processing the signals collected by each buoy 
alone» one can obtain the following noisy measurements: 
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Bearing Indication 



a * 



Bearing is given directly by the buoys or 
preprocessed (filtered) to reduce noise effects. 

Since the bearing is obtained from a vectorial 
sum» it can be reasonably accurate only when the acoustic 
noise is approximately omnidirectional or much weaker than 
the target signal. For some sensors this vectorial sum can 
be limited to a selectable frequency band/ thereby easing 
the noise problem. 

The relationships between a bearing indication 
obtained by a buoy i and the states of the target is shown 
in Fig 3 and given by 

b.Ck) = arctanQy(k) - y.(k)l/(x(k) - x (k)Q t v(k) 



where v(k) must account for all the noise aisturbing this 
measurement/ including acoustic noise/ transducer inaccura- 
cies/ preprocessing noise and errors due to the finite speed 
of the sound in the water. This last factor is caused by the 
time delay between the emission of a sound by the target and 
reception at the hydrophone. At the end of this period the 
target and the hydrophone have slightly different relative 
positions. 



According to 12) tyoical accuracies of _+ 5 
degrees are common for strong sianals and inexpensive OIFAR 
buoys/ and as cited by (31 inaccuracies down to + 1 degree 
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Figure 3 ; Bearing measurement. 
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or less can be obtained with arrays of hydrophones and 
proper preprocessing. 

b. Frequency Indication 

By analyzing records of length T seconds of the 
acoustic signals picked up by the sonobuoys with Fast 
Fourier Transform algorithms^ one can obtain approximate 
frequency spectrum descriptions with a resolution of 1/T 
Hertz. 



By detecting^ measuring and possibly tracking 
some of the strongest frequencies contained in the sig- 
nals very useful data about the target parameters is 
obtained since the received frequencies are functions of the 
originally emitted frequencies^ the soeed^ the heading and 
the position of the targets as shown in the relationship 

f (k) = f (k).c / (c - s (k)) + v(k) (2.25) 

i o i 

where 

s.(k) = -s(k) .cos (b(k) - b (k)l (2.26) 

1 i 

is the relative velocity toward buoy i» as-'shown in Fig 3. 

The term v(k) now has to account also for the 
errors introduced in the various paths of the signal from 
the target to the ship or aircraft# and for variations in 
t arget -sensor geometries during the record time T. 
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As suqgested in and 13] » inaccuracies 

ranging between 0.5 and 0.01 Hertz can be obtained in most 

situations for values of f of the order of hundreds of 

o 

Hertz. 



c. Signal Strength Indication 

The strength of the signals varies not only with 
the strength of the originally emitted sound but also with 
the distance between the target and the buoy. The signal 
strength also depends on the random aspects and the direc- 
tivities of the targets the hydrophones^ and of the 
transmitting and receiving VHF antennas. Random absorotion 
and scattering in the various media from the target to the 
processing equipments may also have a significant influence 
on the signal presented to the processors. 

A mathematical formulation of all these effects 
would be extremely difficult task and would obscure the 
desired distance information. For this reason the informa- 
tion contained in the strength of the collected signals^ 
which can be used by an experienced and skilled audio opera- 
tor» was not incorporated in this study. 

Attenuation factors are thrown away by setting 
the strength of the received signals at a good working 
level . 
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2 » Two or More Buoys 



Considering the signals received by two or more 
buoysf indications can be obtained which may be used in 
addition to# or instead off the measurements available from 
a single buoy. 

a. Frequency Difference Indication 

To process frequency measurements obtained from 
one buoy requires inclusion of the rest frequency as one of 
the states. This is the case because this measurement is 
very sensitive to errors in the assumed rest frequencyf 
i .e . f 



f. = f.C/(C”S.) 

X o 1 



Af. = Af .c / (c - s.) = Af 
1 o X o 

Having the frequencies receivea by two buoys and 
taking their difference one gets 

f. - f. = f .c/(c - s.) - f .c/(c - s.) = 

X j o X o j 

= f .c.(s. - s.) / ((-c - s.).(c - s.)l = 

o X J 1 J 

= f . (s . - s . )/c (2.27) 

o 1 J 



The sensitivity of this difference to errors in 
the assumption of the rest frequency is given by 

A(f. - f.) = Af .(s. - s )/c 
1 J o X j 

3A 




M 



or is reduced by a very large factor since c is of the order 

of 1500 m/s and (s. - s.) is normally very small. 

1 J 

This way» if one has a reasonable approximation 
for f ^ f instead of processing two measurements in a five- 
state plant one can process only one frequency-difference 
measurement in a four-state olant with greatly reduced com- 
puting time and hopefully not a detectable reduction in 
accuracy. 



The variance of the errors in each measurement 
must be added to give the variance of the error in the fre- 
quency difference measurement. 

b. Frequency Ratio Indication 



If one divides the frequencies measured at two 
sonobuoysr a new relationship is obtained that retains the 
information on position/ heading and speed of the target/ 
but , i s independent of the rest frequency emitted 



f. / f.= [f .c/(c - s.)l.C(c - s.)/(f .c)] = 

1 J o 1 JO 



= (c - s ) / (c - s ) 
j i 



(2.28) 



or 



f./ f.= 1 4- (s. - s.)/c . tl + s./c + (s /c)“^ + ...} 

1 J 1 J 11 



= 1 + ( s . - s ) / c 
1 j 



(2.29) 
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As with the frequency difference indication^ 
this measurement can be used in a four-state plant instead 
of using two frequency measurements in a five-state plant. 

If this ratio is generated by the division of 
the frequencies obtained from the analysis of each of the 
signals^ the characteristics of the frequency-ratio measure- 
ment noise is complicated and state dependent. If however 
this relation is obtained as described in the following dis- 
cussion^ this uncomfortable situation can be avoided. 

c. Time Delay Indication 

As shown at the beginning of this section on 
measurement s f the signal collected by a sonobuoy is a noisy^ 
delayedf attenuated and compressed/exoanded reproduction of 
the original signal emitted by the target. The delay is due 
to the finite time the signal takes to reach the buoy and 
the compression/expansion is caused by the variation of this 
delay with time (doopler effect). 

The signals received by two sonobuoys have dif- 
ferent noise contributions and different delay and 
compression/expansion factors even after their strengths are 
equa 1 i zed . 



Supoose one recorded the signals received oy 

buoys i and j from t=ttot=t +T. The signals are 

o 0 

related to the original signal as shown in Equation 
i .e . r 



3b 



r.(t + t') = - d.(t )/(c 

1 • o 1 o o 1 o 



+ c.t'/(c - s.(t ))] 
1 o 



r.(t t') = a.(t + - d.(t )/(c 

J O JO o j o 



+ c.t'/(c “ s.(t ))] 
J o 



s ( t ) ) + 
1 o 



(2.30) 



s . ( t ) ) + 
J o 



(2.31 ) 



If one now amolifies these two signals to a good 
working level/ correlate them by evaluating 



p ( t) = 



r.(t + t').r.(t + t' + t) .dt ' 

1 o JO 



(2.32) 



and look for the value of x that maximizes p (x) one finds 
that : 



(1) if s.(t ) = s.(t ) and/ consequently/ the fre- 
1 o JO 

quency shifts are about the same in both signals/ then r (t 

i o 

t t*) will be simoly an aoproximate delayed or advanced 



replica of r (t + t')/ as can be seen from the above equa- 
J o 

t i ons . 



In this case the maximum of p (x) is easy to 



find and corresponds to the value x at which 

max 



t - d.(t )/(c - s. (t )) + c.(t' +x )/(c - s (t )) = 

o 1 o X o max i o 



= t - d (t )/(c - s (t )) ♦■ c.t'/(c - s (t )) 
o J o jo jo 



solving for x yields 
max 



= td .(t ) - d (t )] / c (2.33) 

max 1 o jo 
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(2) if s.(t ) 4 s.(t )/ d(t ) does not have a 
1 o JO ‘ 

s i dip 1 V “det erm i ned maximum value because the two signals have 
frequency components whose phases are shifted by different 
amount s . 

If ones now change the time scale of one of 

the signals/ say r,(t + t*)/ one has 

J o 

r'.(t’) = r.(t + m.t') = 

J • - JO 

= a. .r fit - d ,(t )/(c - s . (t ) )] + 

J O JO JO 

+ c.m,t'/(c - s.(t ) n (2.3^) 

J o 

Comparing this equation with Equation (2.30) 
one sees that if m can be found such that 

m.c/(c - s.(t )) = c/(c - s (t )) 

JO i o 

o r 

m = Ic - s . (t )) / Ic - s (t )] (2.35) 

JO 1 o 



one again has case (1) and the correlation function will 
give a maximum at approximately 

T = Id ,(t ) - d .(t )) / c (2.36) 

max 1 o JO 

It is interesting to note that Equation (2.35) 
is the same as Equation (2.28). The compression/exoansion 
factor is thus the ratio between the frequencies received by 
buoys i and j . 
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The p ( X » rr ) function/ if the two recorded sig- 
nals were broadband# of long duration and with high S/N 
ratios# would appear as shown in Fig U. 

Low signal to noise ratios could hide the real 
maximum and create false ones# narrowband Or coherent sig- 
nals would present ambiguities with the creation of many 
c 1 ose maxima. 

By considering the availability of eguipments 
and technigues to obtain the proper elimination of the 
difference in freguencv shifts in the recorded signals # 
this study assumed the possibility of having time delay 
measurements available to provide information about the tar- 
get position according to the relationship 

T.. (k) = td.(k) - d.(k)l / c + v(k) (2.37) 
iJ 1 J 

where v(k) must account for the errors caused by the dif- 
ferent noise contributions in each signal# for the accuracy 
and orecision of the cross-correlation algorithm or device# 
and for variations in target-sensors geometry during the 
recording of the signals. 

Inaccuracies anywhere from 0.5 sec down to a few 
milliseconds seem very reasonable to expect in many cases. 

d. Signal Strength Difference 

For the same reasons already explained in 
II#D#l#c the information contained in the difference of 
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Figure 4 



Correlation function. 
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strength of the collected signalSf while of great use for an 
experienced and skilled audio operator/ was not incorporated 
in this study . 



«1 



Ill 



THEORETICAL SOLUTION 



A. RECURSIVE BAYES FORMULATION 

The maximum amount of information at time t, / in a pro- 

k 

babilistic senses is concentrated in the conditional density 
function of the states given all the observations up to this 
time/ o(_£(k)IZ )/ where Z represents the set {z(k)/ z(k-l), 
. . . / ^( 0 ) } . 

By specifying a cost function one determines how this 
information is used to obtain the best estimate [51. For 
example/ the minimum variance estimate is the mean of the 
conditional density. 

From the assumptions made in Section II/A/ the random 
processes x(k) and {x(k)/z(k)> are Markov and the condi- 
tional density can be written in recursive form using Bayes* 
Law/ as is shown in (51 / [61 / [7] / (81 . 

p(x(k)|Z^) = c(k) .p(x(k) I.Z’^"^) .p(z(k) ! x(k) ) (3.1) 



where 



p ( X ( k ) { Z^~h 



p(^(k-l ) I Z^"^) .p(x_(k) !^(k-l ) ) .d^fk-l ) 

(3.1a) 



and 



1 / C ( k ) 



= p(2(k) = 



p(x(k)!Z^ ^.p(z(k) 



(3.1b) 



Given the initial condition p(j^(0)|Z~^) = p(x(0)) and the 
probability law of the system random input and measurement 
noise/ p(v^(k)) and p(^(k))/one can proceed to evaluate these 
equations and the fflterinq problem can be regarded as hav- 
i ng been solved when p(x^(k)!Z ) can be determined for all k. 

For any situation different from the one where the sys- 
tem and measurement equations are linear and the apriori 
distributions are Gaussian, the density p(x^(k)lZ ) is not 
Gaussian and generally cannot be determined in closed form. 

In the problem studied here/ the conditional densities 
in Equation (3.1) are nonlinear functions of the states and 
the measurements. Linearizations/ Taylor expansions and 
Gaussian approximations will be used, however, as an 
engineering comoromise, in order to avoid complex numerical 
integrations in n dimensions. This aoproach was chosen tak- 
ing into account the limitations on processing eouioment 
available for the practical solution of the addressed track- 
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ing problem 



B. ONE-STEP PREDICTION 



Let us assume that at time ' during a filtering pro- 

cess/ an approximately Gaussian density p ( x ( k- 1 ) ' Z has 
been obtained with mean value p (k-1) and covariance matrix 
M(k-l) . 

If one chooses the mean of this conditional density as 
the estimate at time ^/ 

x_(k-n = E tx_(k-l) =_p(k-l) 

the estimation error at time k-l has the same probability 
density as p(j^(k-l)lZ^ ^ but with zero mean/ i.e./ 

E[e(k-1)) = El(£(k-1) - ^(k-l))|Z^~^ = 0 (3.2) 

E(e(k-l)e'^(k-l)] = P(k-l) = ►^(k-1) (3.3) 

At time t, a new set of measurements arrives and a pred- 
k 

iction at that time must first be obtained based on the 

measurements up to t . 

k-l 

The dynamics of the system are described by Equation 
(2.3) and it can be seen that/ even if p ( x ( k- 1 ) | Z^~^) were 
really Gaussian/ the density p (_^( k ) | Z^~^) would not normally 
be because of the non 1 i nea r i t i es involved. Nevertheless it 
will be assumed that a reasonably close approximation to a 
Gaussian density exists and proceed to find approximate 
values for the first two moments of the conditional one-step 
pre-diction density p ( x ( k ) ! z'*’^) . 



The first step is to exoand Equation (2,3) in a Taylor 

series about the estimate x(k-l). The dependence on t, and t. 

— k k- 1 

\ 

k— 1 

is dropped for simplicity. Note that Z i-s known and was 
used in determining ^(k-1). 

x(k) = f(x(k-l)) + g ( X ( k - 1 ) ) . w ( k* 1 ) = 



f(J(k-U) + 

— — 3x 



. tx(k-l ) 
x(k-l) 



- x(k-l)J + 



+ q(x (k-l ) ) ,w(k-l ) + 



_f (£( k- 1 ) ) - ^ 

i=l 



n n 

^ I I 

i=l j=l 



3f_ 

X 



• e , ( k- 1 ) 
x(k-l) 



~ 




3f 


,e.(k-l).e. (k-l) 

1 j 


3x. 3x. 
1 J 

w 


x(k-l^) 



f q(x(k-l)).w(k-l) - 



n 

I 

i=l 



"'Si 






, w ( k- 1 ) 


, e ( k - 1 ) + 


dx 


i 


— 


- 


i(k-i) J 




individual components of x 



(3. a) 



g. is a column vector of g ! and w is a component of w. 

X — i “ 

1 . Linear Approximation 

In this case only the terms in the expansion which 
are linear in e(k-l) and w(k-l) are considered^ i.e.^ 

k ) = _^(^(k-l)) - ($1 ( K- 1 ) . e ( k- 1 ) + 

+ q(x(k-l)).w(k-l) 



^5 



(3.5) 



where 



<J)(k-n = -^ f_(i(k-n (3.6) 

From Equation (3.5) the mean value of p ( x ( k ) ! )» 

which is chosen as our predicted value» is 

i(k}k-l) = E(_x(k)!Z^"^ = Xflfk-D) - 

- <j> (k-1 ) .E te(k-l ) } + g(^(k-l ) ) .E (w(k-l ) J 

But e_(k-l) and ^(k-1) have zero means» therefore^ 

J(k!k-1) = f(i(k-l),t, ,t, J (3.7) 

~ — — k k-1 

Thus the one-step prediction error has the mean 
E[e(k!k-1)) = E[(£(k|k-1) - ^ ^ = 

= E [(j>( k-1 ) .e(k-l ) - g(i (k-1 ) ) .w(k-l )J = 0 

The prediction error covariance matrix, which is 
also the covariance matrix of p(x(k)| Z^~ ^), is, using simpli- 
fied not at i on , / 

P(klk-l) = E(e(kS*k-l)e^’^(k!k-l)l = 

= E (k-1 ) . 

TT TT T*T* TT* 

= E l ^ . (j) - (j>.e._w .q - g.w.e.(j) g.w.w'-.g-'-] 



e(k-l) - g(x(k-l)).w(k-D) I 



...]TJ = 



Ub 



From the assumptions about the independence and 
discrete white noise characteristics of w, the terms 

(j) (k-l).E[e(k-l).w'^(k-l)] = 0 

j(i(k-l)).E[w(k-l).e^(k-l)) ,<j)'^(k-l) = 0 

and 

P(k|k-1) = (j) (k-l).P(k-l). (j)'^(k-l) + 

+ g(J^(k-l ) ) ,Q(k-l ) .g^ (£(k-l ) ) (3.8) 

2 . First Order Approximation 



In this development one retains all terms in Equa- 
tion (3,4) up to first order partial derivatives^ i.e.» 



_x ( k ) = ^(j^(k-l)) - (j) ( k - 1 ) . e^( k - 1 ) + 

+ £(j^( k- 1 ) ) ,^( k- 1 ) - 

m 

- y A,(k-l).w (k-1) . e(k-l) (3.9) 

i=l ^ 

where 

A, (k-1) =-^ a. (J(k-l),t ,t ) (3.10) 

3x — k k-1 



As in the last 
^(k|k-l) = 
+• g (_x 



section/ the predicted value is 

l(^(k-D) - <j) (k-1) .E le(k-l )] + 

m 



(k-1) ) .£ (w(k-l )] 



- I 

i=l 



A .EIw (k-1 ) .e(k-l )] 
11 — 
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From the assumptions about ^(k-1) and e(k-l)» 
predi ct i on is 

_x(k|k-l) = ^(^(k-1)), 

the same as in Equation (3.7), 



The one-step prediction error has zero mean 
covariance matrix given by t in simplified notation# 



[' 



P(klk-l) = E (<j,(k-l).e(k-l) - g ( x ( k - 1) ) . w ( k- 1 ) + 



m 

+ y A. (k-l ) .w . (k-l ) ,e(k-l )1 [ 

1 — 

i=l 



= E [ <j), e_.£^ . (j) - (j) . e_, ^ (j) . e . e'*' . A-^ . w _ - 



T „T 



m 






.T aT 



i=l 



1 1 



T T 



+ y A.,e.e .d) .w. - y A ,e,w ,q .w + 

i=l " - - ^ 1 i=l " - - i 

m m 

+ y y A,.e.e^ .A*y.w..w ) 
i=ij=i ^ ^ j 



m 



” d_»_w • £* • (j) " g.w.w'T.gT - g.w.^.A^, 

m _ _ m 

.T 



From the assumptions about w# one has that 



E (e(k-l ) .w^(k-l)) = 0 



E Ie(k-1 ) .e^ (k-l ) ,w^ (k-l)l = E te ( k- 1 ) . eM k- 1 ) 1 . 



E (w (k-l )) = 0 

i 



E[^(k-l).e (k-l).w (k-l)) = E(w(k-l),w (k-l)l. 

~ i — i 



E (eMk-1 )] = 0 



the 



and 



w + 
i 



as 



t[e(k-l).e ( k - 1 ) , w . ( k - 1 ) , w . ( k- 1 ) ] = 

- - 1 j 



= EIe(k-l).e^(k«l)] .Elw. (k-l).w. (k-1)] = 

- - 1 j 



= P(k-1 ) , q. . (k-1 ) 



where q..(k-l) is the jflth term of Q(k-l) 



Thus f 



P(kJk-l) = (j) (k-1 ) .P(k-1 ) .(() ^ (k-1 ) t 



+ 2(^(k-l ) ) .Q(k-1 ) .jMjc (k-1 ) ) t 

T 



m m 

+ y y A .P(k-1 ) .At .q. . (k-1 ) 



=1 j=l 



(3.11) 



If the Q matrix is diaqonal# q.. = 0 for j t 1 and 



P(klk-l) = <<)(k-l).P(k-l).4,^(k-l) + 



+ g(x(k-l)).Q(k-l).qMx(k-l)) + 



m 



y A. .P(k-1).A^ .q. (k-1) 
X X xi 

x=l 



(3. 12) 



3 . Other A oproximations 



Depending on how many terms of Equation (3.^<) are 
considered^ other different prediction values and prediction 
error covariance matrices are obtained. Higher moments of 
the conditional densities would be required with increased 



complexity and computational burden 



C. UP-UATING THE ESTIMATE 

Let us assume that we have an approximately Gaussian 
k-1 

density p(x.(k)lZ ) whose mean value is the prediction 
x_(k|k-l) and whose second central moment is the prediction 
error covariance matrix P(k!k-1). 

A new set of measurements is given# consisting of non- 
linear noisy observations of the states# represented by the 
re 1 at i onsh i p 

k ) = ^ ^ ^ k ) 

It is now necessary to process z(k) and extract the 
information it brings. A new approximate Gaussian density 
p(j<(k)}Z ) is assumed to be generated whose moments are the 
new estimate £(k) and the estimation error covariance matrix 
P(kl . 



This new density is given by Equation (3.1) or by 
p(£(k)JZ^) = p(£(k) |x_(k) ) .p(_x(k) :Z^~1) / p(£(k) ! 

(3.13) 



p(_x(k)5Z^^ is assumed to be approximately Gaussian with 
moments ^(kJk-1) and P(k|k-1). 



Expanding the measurement eauation around x(k|k-l) one 



has 



i(k) = h(x(k;k-l)) + 






.(x(k) - _x(kjk-l)] + ... + ^(k) 

x(k|k-l) 
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Retaining only the linear terms^ one has 



^(k) = ^(jT(k!k-l) + H(k).[j^(k) - x(k!k-l)] +_v(k) 

(3. U) 

where 

H(k) = ^ h(i(k ! k-1 ) , t, ) (3.15) 

9x — — k 



From Equation (3. 


,1^) one obtains 


E lz_(k) ! 


= h(£(k|k-l)) 


and it is easily 


shown that 


Var (z(k) !Z^~ 


h = H(k) ,P(k ! k-1 ) .hF (k) + R(k) 



Also from Equation (3,1^)/ 

E (z^(k) ;^(k)J = Jh(^(k{ k-1 ) ) f H ( k ) . [_x ( k ) - £(k!k-l)] 

Var(^(k)!j<(k)J = R(k) 

If one now assumes that p(z(k)'^~^) and p(z(k)lx(k)) are 
approximately Gaussian densities with the moments given 
above/ then Equation (3,13) becomes 



p(£(k)',Z^) : 


T, 1/2 

: |H(k) ,P(k ! k-1 ) ,HMk) + R(k)| / 


/ 


^/9 1/2 1/2 
((2Trr^.|P(k!k-l) 1 ,|R(k)| ] . explB) 


where 
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explB) = exp 



{■ ’ ['* 



(x (k) - i (k ! k-1 ,P ^(k 1 k-1 ) 



. ( . . . ] J 



t_ 2 (k) - ]i(£(k;k-l)) - H(k).[_x(k) - i(k!k-l)J]T. 

R" 



■ [. 



■^(k) . 1. . .]J - 



rz(k) - ^(^(k|k-l))]T. 

tH(k).P(k|k-l).HT(k) + R(k)f^.t...] 

(3.16) 



Simplifying the notation^ this exponent can be rear- 



ranged in 


the 


f 0 1 lowing 


way 




T 

( X - X ) . 


P~\ ( X 


- 5) + (2 


” h ( X ) 


- H. (x - J .R"\ 




- (z 


- h(x)J^. 




■ RP^ (z - h (;)] = 


= (x 


r 

- x) 


^P~\x - x) 


+ (x - 


i)'^H'*^R"^(x - x) - 



- l2 - h()?)f FT^Hfx - x) - (x - i)'^H^R'^(z - h(i)l + 

+ ( z - h ( i R“^[ z - h ( X ) I - 

- tz - h(i))^ [HPH^ + Rf^(z - h(i)I = 

^ T -1 T -L 

= (x - x) (P + H R>l)(x - x) - 
^ T T -1 

- (x - x) HR Iz - h(x)] - 

- (z - h(x"))^R“^(x - x) - 

- (2 - h(^)l^CR“^- (HPh'^ + R)“^] [z - h(x)J 
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Using the matrix inversion lemma (91, 

- (HPh'^ +R)"^= R H(P"^+ 

Since R and P are symmetric matrices, 

R~^(P~^+ H’^R"^fVR~^= A^(P~^ + hfp'^lA 

where 

A = (P"^ + h'^R^Hf^^R"^ 

and the exponent of Eauation (3.16) becomes 
( X - i )^ ( P“^ + R"^H ) ( X - i) - 

• (x “ x)^hFr^ 2 “ h(x)l • 

- (z - h(x)] R“Vl(x - x) - 

- ((P"^+ z - h(i))J^(P“^t h^fT^ll...] 



which is a quadratic form that can be expressed as 




(P“^ + h'^ R“^r^^R"^lz 



h(x)] 



3 - 



(P“^ + 



R"^H ] 



And finally the conditional density of Equation (3.16) 
can be put in the form 

o(x(k)!Z^) =a(k).exo |^-^(ii(l<) - ^ ( k ) 1^. P"^( k ).(... 1 
where, by definition. 
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x(k) 



= X ( k I k “ I ) 



P(k) .HT(k) .R-^k) . [z(k) 



h(x(k|k-l),t )1 
— ~ k 



(3.17) 



and 

P"\ k ) = P"^ k S k - 1 ) + ( k ) R”\ k ) H ( k ) 

OPf using the matrix inversion lemma again^ 

P(k) = P(klk-l) - P(k|k-l)H^(k). [H(k)P(k!k-l)H'’^ (k) + 

+ R(k)]"^H(k)P(k I k-1 ) 



or 

P(k) = Cl - G(k)H(k)J .P(klk-l) (3.18) 

where 

G(k) = P(k|k-n.H^(k). [H(k).P(k!k-l).H^(k) + R(k)J 

(3.19) 



From hquation (3.19) one obtains - 
G(k) = P(k!k-l).HT(k).R"^k). 

(H(k) .P(k I k-1 ) .H^(k) ,R"^k) + 



Premultiplying by P(k)P^(k)^ and using the expression found 
earlier for P \k) gives 



G(k) = P(k).(I + H^R~^HP(k ;k-l )) H^R"\ 

(HP(k 1 k-l)H^R"^ + I)"^ = 



5a 



= P(k).tH^R"^+ H^R^HP(k ! k-1 )H^R ^ . 

(I + HP(k| k-1 )h’’^R'^~^ 



and f i na 1 1 y » 

G(k) = P(k) .H^(k) .R"\k) (3.20) 

Applying this result to Equation (3,17) gives 

5_(k) = _x(k|k-l) ♦ G(k).[j^(k) ~ Ji(_x(k!k*l)>tj^)) 

(3.21) 

If more terms were taken in the expansion of h(x(k)»t, ) 

— — k 

than the ones considered in Equation (3.1^)» different esti- 
mates and covariance matrices would be obtained. Higher 
moments of the conditional densities would then be required. 

iNe are now in position to restart the one-step predic- 
tion of Section III^B and process a measurement ocurring at 

The initial step in the whole process is at t where it 

o 

is assumed that an approximate Gaussian distribution exists 
for the initial value of the states. j^(0). Choosing the ini- 
tial estimate as the mean value of this distribution, one 
has 



J(0) = Etx(O)) 



Ele(0)l = E(i(0) - _x(0)l = 0 
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P(0) = E (e(O)e^^(O)] = Var[^(0)J 

The set of Equations ( 3 , 7 ), ( 5 . 8 ), ( 5 . 18 ), (3,19) and 

is generally known as the Extended Kalmam Filter 
equat ions. 

D. PARAMETERS FOR THE TRACKING PROBLEM 

In order to apply the above relations to the passive 
tracking problem^ a series of vectors and matrices first 
must be determined, 

1. System Dynamics 



Applying definition (3.6) to Equation (2.10) yields 



(J)(k) 



(3.22) 



1 



0 



At .cos b ( k ) 



At,s(k),sin b(k) 



0 



At . s i n b ( k ) 



At ,s(k) .cos b(k) 



0 



0 



1 



0 



0 



0 



0 



1 



where At 




Breaking up Equation (2.11) into its column vectors 



gi ves 
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g (x(k)/t »t ) = 
-1 -• k+1 k 



— At^ .COS b ( k ) 
2 



— At^.sin b(k) 
2 



At 



(3.^3) 



g (x(k)»t, rt) = 
_2 — 1 ^+]^ k 



•^At^.s(k),sin b(k) 
2 



■^At^.s(k).cos b(k) 
2 



(3.2a) 



At 



Definition (3.10) requires 



(k) 



(3.25) 



0 0 0 
0 0 0 
0 0 0 
0 0 0 



— At^ .sin b ( k ) 



-i-At^ . cos b ( k ) 
2 
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A2(W) = 

0 0 -^At^.sinb(k) 

0 0 .cosb(k) 

0 0 0 

0 0 0 



(3.^6) 

- — A . s ( k ) .cos b(k) 

2 

At^.s(k).sin b(k) 

2 

0 

0 



Assuming independent random forcing inputs yields 



Q(k) = E(w(k)w (k)) = 



onk) 0 
s 



0 o;^(k) 

b 



(3.27) 



Equations (3.8) and (3.11) require the product 



g(x(k)»t ft ) . Q ( k ) . g^ ( X ( k ) f t ft ) = 
k+1 k - - k+1 k 



= At^ . 



o (><) 0 

s 



o^(k) 

D 



(3.28) 



where 



a = 



. (cos^ b(k) . a^(k) + s^ ( k ) . s i n^b ( k ) .^ ^( k ) ] 

S D 



b = -rAt^.sin b(k).cos b(k).(o^(k) - s^(k).Q^(k)) 
4 s ^ b 
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I# 



c = yAt Is i n^b ( k ) . a^( k ) + s ^ ( k ) .cos^b ( k ) . ( k ) ] 

H s b 

and 

d = -rAt.cos b(k),a^(k) 

2 s 

e = - 4 At . s ( k ) . s i n b(k).a^(k) 

^ b 

i = -^At.sin b(k).a^(k) 

2 s 

j = ^At.s(k).cos b(k).a^(k) 

Z b 



If the expanded state variable formulation is 
the matrices given below are needed 



0®(k) = E [w®(k)w®^(k)l = 



Q(k) 



I 

I 

I 

I 

I 

I 0 

I 



(3.29) 



® n eT _ 
g .Q .g = 



g.Q.g 



.2 2 

A^ • 



(3.30) 






0 ^ 1 



(3.31 ) 



used / 
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I 



I 



2 . Measurement s 



In the following development Equation (3.15) is 

/ 

applied to each of the measurements described in Section 
II/D to yield the appropriate linearized observation matrix. 

a. Bearing Measurement 



From Equation (2.24) one has 



3 b. = (x - x.)^/((x - X.;? + (y - y.)^l . 

55 111 1 

-^((y - y. )/(x - X.)) 

95 1 1 



and then 



1,1 



(k) = 



b^(^(k|k-l) ,t^ ) 






ly(kik-l) “ y . ( k ) ] /a tx(klk-l) 

1 



- x^ (k)l / aj 



0 0 
(3.32) 



where 

a= d?(k) = (x(klk-l) - x,(k))^ + Iy(k!k-l) - y (k)]^ 

11 i 

(3.33) 

and d^(k) is the predicted distance of the target from buoy 

i . 



b. Frequency Measurement 



From Equation (2.25)» 
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i 




(c/ (c 



f . 
1 



s.)].Af + (f ,c/(c - s.)^].As. 

1 ag o o 1 9^1 



f /f 

i o 



J- f 
3C o 



/ ( f .c) .As 
1 o 9C i 



where s is given in Equation (2.26) and 
i 



s . = - cosCb - b.).As + 
9C 1 1 9C 



+ s.sinCb - b ) . -^ b - 
i 9C 



- s.sin(b - b ). — b 
i 95 i 



and then/ 



H^.(k) = (x(k!k-l)/t, ) 

fi 9x 1 — k. 



(3.3a) 



■t 



ly(k}k“l)-y^(k)]/5 -(x(k!k-l)-x^(k)J/^ 3 Y 



f /f 
i 



where 



(3.35) 



f = f . (x(k ! k-1 ) , t ) = f (kik-n.c / [c + 
1 1 — k o 



+ s(k!k“l).cos(b(ki'k“l) • b (k))) 

i 



[' 



b.(k) = arctan [y (k I k-1 ) - y (k)]/ 
1 ■ 1 



(x(k!k-l) - X (k)J 

i 



(3.36) 



Y- f ^ . s ( k 1 k- 1 ) . s i n (b ( k 1 k- 1 ) - b (k)) / 
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(f (k!k-l).c) 
o 



e = - ff. cos(b(k|k-l) - b. (k)) / If (kl k-l).c] 



6 ” a^Y ' o (3«33). 



c. Frequency Difference Measurement 



From Equation (2.26)^ 



— (f. - f.) = f /c . — { s . - s . ) 

95 1 J 0 35 1 J 



where — s is given in Eauation (3,34). Thenr 
95 1 



H ^.(k) = A (f - f ) 
di • 3x1 j 



where 



= f /c . (h h h h ] (3.37) 

o 12 3 4 

x(k|k-l) 



h = a .(y(klk-l) - y (k)3 - a .ty(k'k-l) - y ( ic ) ) 

1 i i j j 

h = a .(x(klk-l) - X (k)3 - a .(x(k!k-l) - x (k)3 

2 J j i i 

h = cos(b(k|k-l) ~ b (k)) - cos(b(klk-l) - b (k)) 
2 j i 

I 

h, = s ( k ! k - 1 ) , I s i n (b ( k ! k- 1 ) - b (k)) - 

sin(b(k'.k-l) - b (k) )3 

j 

a. = s(k{k-l).sin(b(kjk-l) - 6 (k)) /„ 



b.(k) as in (5,36) and n as in (3,33), 
1 
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d 



Frequency Ratio Measurement 



From Equation (2.28)/ 



-L (f./f, ) = 1/c 
yC • 1 j 




s ) 
j 



Comparing with the frequency difference measure- 
ment / one has 



H (k) 
-ri 



6 • 



= (f /f ) 

3x i j 



= 1/f 



x(k|k-l) 



Time Delay Measurement 



H (k) 
dl 



(3.38) 



From Equation (2.36)/ 

-i- T = 1 /c . -i-(d - d ) 

ij 3^ i j 

where 

d^ = (x - X )^ + (y - y )^ 

i i i 



and then/ 



H,..(k) = -It.. = 1/c . Ih , h 0 0) (3.39) 

•TiJ 3x iJ ^ 12 

x(k|k-l) 

where 



h,= (i(Wlk-l) - x.(k)l / d.(k) - 
1 11 

lx(k!k-l) - X .(k)] / d .(k) 

J J 
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h = [y (k ! k-1 ) - y (k)J / d (k) - 
^ i i 

ty(klk-l) - y.(k)J / d.(k) 

J 

and d is given in Equation (3.33). 
i 



bU 



IV. SPECIAL CHARACTERISTICS OF KALMAN FILTERS 



Equations (3.7) and (3.8) characterize the one-step 
prediction of an Extended Kalman Filter; Equations (3.18)» 
(3.19) and (3.21) represent the estimation update. 

When implemented in a computing device# the timing of 
these steps is shown in Fig. 5. The actual state is 
represented by a broken line and the output of the filter by 
the solid line. In general# one is dealing with vector 
processes so Figure 5 is r ep r esen t a t i ve of only one com- 
ponent. Nevertheless# this, reoresentation is heloful in 
visualizing the steps involved in the estimation process. 

In a general problem the times of occurrence of meas- 
urements are not equally separated and are not known in 

advance. At a time t # when a new set of a measurements 

k-1 

becomes available# the output of the filter is still the 

last estimate made close to t # unless predictions are made 

k-2 

at regular time intervals independent of the spacing between 
measurement s . 

A time period T is spent in a prediction phase to 

P 

evaluate Equations (3.7) and (3.8)# after which the best 
information about the state of the plant is the oredicted 
value x(k-lIk-2). 
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Figure 5 ; Timing in a Kalman Filtering Process. 



If the times of occurrence of measurements were known 



in advance^ the predictions 
tern was idle waiting for the 
be made zero. 



could be made any time the sys- 

new measurements and T could 

P 



Nextf the q measurements are processed using Equations 

(3.16)f (3,19) and (3.21) and a time T is spent which is 

e 

relatively large mainly because of a matrix inversion of 

dimension (q x q)^ in the gain equation^ (3,1,9), Only after 

this time has passed is the new estimate x(k-l) available. 

This estimate will be the output of the filter until t, when 

k 

the process is repeated. 

Some important points in this estimation process must 
be stressed: 

-** the new estimates are only available T + T 

P e 

seconds after the arrival of new measurements. 

--- a major oart of the time period of duration T is 

e 

spent in the gain equation due to the required matrix inver- 
sion, 

--- from the enq of the computation of an estimate 

until the arrival of a new set of measurements the computing 

equ ipment is idle and thus free to execute other chores. 

--- the one-step prediction is based on a 1 i near i zat i on 

of the plant dynamics over an extended periods for example^ 

from t to t , 
k-1 k 

--- the measurement equations are based on a lineariza- 
tion about the predicted values. 
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i. 




A. PARTITIONING OF MEASUREMENTS 

When some of the measurements come from statistically 
independent sources/ steps can be taken to alleviate orob- 
lems which arise. Considering the special case when all the 

measurement noise components are independent/ the covariance 

\ 

matrix R(k) is diagonal and the following development 
app lies. 

1 . The Linear Case 

Assume a dynamic system with linear observations of 

the form 

^(k) = ,H(k) .^(k) + V (k) (a. 1 ) 

where H(kJ is not a function of the states. 

The estimation equations of Chapter III become the 
standard Kalman Filter equations 

G(k) = P{k|k-l)H^(k) [H(k)P(k!k-l)HT(k) + R ( k ) (a. 2) 

^(k) = i(k!k-l) + G(k)(_z(k) - H(k).^(k!k-l)l (^.3) 

« 

P(k)=lI-G(k)H(x)]P(k;k-l) (a.y) 

The measurements are assumed to occur simultaneously 
and/ as shown for the first time in (91/ for this linear 
case/ the results obtained by processing all q measurements 
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at once using these eouations are the same as the final 
results obtained by processing one measurement at a time. 

In the latter approach^ the result of processing one 
measurement component is used in the following computation 
to process the next measurement component. 

A general way to show that the two approaches yield 
the same estimate after all measurement components have been 
processed is given in the following paragraphs. 

In short it is to be proved that the estimates x ( k ) 
and P(k) obtained by the use of Equations (^.2)/ (4.3) and 

(4.4), when all measurements are grouped in a vector z, are 
the same as those obtained by the use of the following 
iterative equations q times, once for each componentJ 

G . = P. ,.H J/ [H. P. hT t R. J (4.5) 

1 x-1 1 1 1-1 1 1 

X. =x +G,(z, “H.x, 1 = 

—1 —1-1 1 —1 1—1-1 

= (I - G,H. ) X . t G z (4.6) 

1 1 —1-1 i i 

P. = [I - G. H .] .P , (4.7) 

1 11 1-1 

i - If 2, .../ q 



where 

X = x(k!k“l) 
— o — 



P = P(k ! k-1 ) , 
o 
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P = P(k) , 
q 



X = X ( k ) 
-q - 





^ z 

1 




r H n 
1 








^2 




H2 




''2 


z ( k ) = 


• 




• 


X ( k ) + 


• 


\ , 


• 




• 




• 




z 




H 




V 




_ q _ 




- 9- 




- q - 



and E tv 1 = 
i 



If Equation (q.6) is 
relationship between ^(k) 
resu 1 t s 



applied a times» the following 

(or X ) and x(k'k-l) (or x ) 
- q - -o 



x(k) =tI-GHltI-G H ]...tl - G H ]x(k!k-l) + 

— q q q-1 q-1 11 — 



+ Gz +tI-GHlG,z +...+ 

q q q q q-l q-l 



+ tl - GqH^tl - . . . 1 1 - G^H^IG, z. 



2 2 11 



(y.8) 



and recursive use of Equation (y.7) gives 

P(k) = tl - G„H^] (I - G ,H J...tl - G,HJP(k|k-l) 

q q q-l q-l 1 1 

(y .9) 



Comparing (y,8) to (y.3) and (y.9) to (y.y) one sees 
that the assertion can be proved by showing the validity of 
the relations 
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II - G(k)H(k)] = II -G H ] [I - G H )..,ll - G H 

q q q-1 q-1 1 1 

(a. 



G(k) 



(n X 1 ) 


1 

1 

1 

1 


1 

1 

1 

I 


(n X 1 ) j; 

1 

I 


II - G H 1 ... 1 1 


1 

- G, j . . . 


1 

1 

1 


1 

(I - G H JG J 


q q 


2 2 1 j 

1 


1 

1 

1 


q q q-i| 
1 


UK 


1 


1 


1 



or# in 3 compressed way^ 



G(k) 



( n X (a- 1 ) ) 



( n X 1 ) 



( I - G H ) G* 

q q 




(a. 12) 



and 



[I - G(k)H(k)] = [I - G H 1(1 - G*H *1 (a. 13) 

q q 



where G is the gain obtained from Equation (3.19) if 
has a q"l dimensional measurement vector# and H and 
equivalent matrices# as shown below 



H(k) 



H 

q 



R(k) 



0 





(a. la) 



(a. 15) 



1 

10 ) 

X 1) 
G 

q 

1 1 ) 



one 
R are 
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G* = P(k!k-1)H*^ (H*P(k I k-1 )H*^ + R* ] ^ (a. 16) 



Let us expand Eauation (U.2) and try to manipulate 
it into the form of Eauations (4,10) - (4.13). We have 



[HPH + RJ = 





R^l H^PH^ 

1 

1 

1 


1 • • • 

1 

1 

1 

• *1 

1 


1 

1 

1 

H 

1 


H PH^ 

1 q 




H2PH^ 

1 


^ 2 \ 

1 


1 

1 

1 

1 

1 


H-PH 
2 q 




1 

1 

1 

1 


1 

1 

1 

1 


1 

1 

1 






• • • 1 • • • 

1 

1 

1 

1 


i • • • 

1 

1 

1 - 
1 


1 

1 

1 

1 


• • • 




H PH?^ ! 
L q 1 1 

1 


1 

1 

1 • • • 

1 

1 


1 

1 

1 

1 

1 


H Plf + 1 

q q 



or 



(HPH" + R] = 



T 

-3 



^2 



( 4 . 17 ) 



where 



h*ph*t ^ R* 



1-2 



H*PH^ 



5-3 - « PH - S.2 



a , = H PH" + R 
4 q q q 



Next define the scalar c as 



c = a 



- a^/Tl. 



3 12 
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1 



Thus < 



C = H Ph'^ + R - H PH*'^ [H*PH*'^+ R*r^H*PH'^ 

q q q q q 



Using the definition of G in Equation (q.l6) gives 



c = H Ph'^ + R - H G*H*Ph'^ - 

q q q q q 



= H [I - g*h*]ph'^ + R 

q q q 



(q.i8) 



It is shown in CIO] that 



A, la" 


“X 


OD 

cr 


1 1-2 
1 

«» W 1 «» 




1 i -2 

«» J w 


1 

T 1 

3 ^ a , 




b ^ i b, 


-3 1 4 




-3 4 






L ~i 



where 



B, = 



-2 









^3 



T.-l -1 



b, = c 



_ -1 



(a. 19a) 
(4.19b) 
(4.19c) 
(4.190) 



Using these re 1 at i onsh i ps > 



G = Ph'^ CHPh'^ + R)~^ = 



= CPH*'^ I PH ] 

q 



^3 



^2 
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= [PH*^B, + PH^b^ I PH*^b + PH^b ] = 

1 q-3 - 2 q 4 

= (C I 0] (a. 20) 



Applying Equations (4,lPa) - (4,19d) to (^.20), 
D = - PH*'^ tH*PH*'^ + R* r^H*PH'^c"^+ Ph'^C~^= 

q q 

= - G*H*PH^c"^t PH^c“^ (I - G*H*]PH^c"^= 

q q q 

= Cl - G*H*JPH^[H (I - G*H*)PH^ t R 

q q q q 



Let P* = II - G* H* ]P be the estimation error 
covariance matrix after the processing of a-1 measurement 
components^ as seen from Equation (3.18). Then> 



0 = P* H^CH P*hT t R = G (4.21) 

q q q q q 



is the gain when processing the qth measurement. 



We also have 



C = PH*^(H*PH*^ + R*r^ + PH*^CH*PH*^ + R*r^H*PHTH PH*T 

q q 

(H*PH*^ + R*r^c~^- PH^H PH*T CH*PH*T ^ R*rV^ = 

q q 



= G* + G*H*PhTh G*c-1- PH^H G*c-1 

q q q q 



Since c is a scalar^ one can write 



C = G* - (I - G*H*1 PH'^c'^H G* = G* - DH G * 

q q q 



and/ therefore/ 
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C = (I - G H JG 

q q 



(a. 22 ) 



and from (^1.15)/ 



= [' 



G = I tl - G H J G 

q q 



,] 



(a. 23) 



It is also the case that 






* . 



GH s fl - G H ] G G 

q q q 



] 



H 



q^ 



= tl - G H ]G*H*+ G H 

q q q q 



so 



I-GH = I- GH-II-GH]G*H* 

q q q q 



andf thus 



I - GH = [I - G H ] [I - G*H*) 

q q 



(4.24) 



Equations (4,23) and (4,24) are the same as Equa- 
tions (4,12) and (4,13)# proving the initial proposition. 

Computationally# one can replace a (qxq) matrix 

\ 

inversion and a (nxq) x (qxq) matrix multiplication by only 
q scalar divisions. For an indication of how much computing 
time can be saved# consider that a matrix multiplication 
requires approximately nq^ scalar multiplications and addi- 
tions# and that a matrix inversion requires at least q mul- 
tiplications and additions# not counting logic and indexing 
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time nor considering storage requirements. 

The partitioning of the measurements then 

2 

corresponds to trading at least q (n + q) multiplications 
and additions for q divisions. For a machine with average 
multiplication time of 10 m i c roseconds » addition time of 2 
microseconds and division time of 12 m i c roseconds f if n and 
q are equa 1 to ^ a saving of about 1.5 ms is obtained; if n 
and q are of the order of 20^ this saving could be of the 
order of 200 ms at each measurement time. 

Besides this comouting time saving^ the processing 
of one measurement at a time orovides partial estimates that 
can be used as improvements over the predicted values* as 
shown in Figure 6. Also* the final estimate is obtained 
sooner* although with the same accuracy as when the measure- 
ments are all processed simultaneously. 

2. The Nonl i near Case 



We assume the same system as before but with non- 
linear observations of the form 

\ 

^(k) = 



By linearizing this equation as shown in Chapter 
III* using the definition 

H(k) = h(x(k|k-l)*t ) * 

3x k 

the proof of the last section is still valid for the 
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q measurements. 



available . 




/x(k) 



x(k) = x(k) 




first measurement 



Figure 6 



Partitioning 
Linear case. 



of measurements 
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linearized equations 



Iff however^ the partial estimates obtained by pro- 
cessing one measurement (hopefully a better approximation to 
x(k) than x(k}k-l) results) are used to calculate the H 
vector for the processing of the n*xt measurements better 
accuracy should be expected than when all H, vectors are 
obtained using the relatively crude predicted value* This 
effect should be particularly significant when the time 
between each group of measurements is relatively large. 

Thus# partitioning the measurements in an Extended 
Kalman Filter should provide improvement in estimation accu- 
racy in addition to the advantages pointed out in the previ- 
ous section. This improvement can be visualized as shown in 
Figure 7, 



The expected improvement in the estimation process# 
brought by partitioning the measurements in Extended Kalman 
Filters# can be explained in terms of the conditional densi- 
ties. 



Consider two independent simultaneous measurements# 



z, (k) 




h.(x(k)) + V (k) 


1 




1 

! 1 


1 






h. (x(k) ) + V (k) 


2 J 




1 

CM 


2 j 



The conditional Density after processing the two 
measurements is# from Equation (3,1)# 
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x(k) 




measurements 
one at a time 



alltogether . 



Figure 7 



Partitioning of measurements - 
Nonlinear case. 
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p(2(k)!Z^) = p(_x(k) (k),z^(k),Z^"^ = 

= p(z^ (k),Z 2 <k) l^(k)).pU(k) :z‘^"^)/p(z^ (k),z^(k) IZ^"^ = 

= [p(z 2 (k) ! z j^(k) ,x(k) ) / p(z^ (k) Iz^ (k) ,Z^"^)1 . 

[p(Zj^ (k) ! x_(k) ) .pU(k) lZ^~^/p(z^(k) |Z^~^)] 

But f f rom (3.1) 

p(^(k) : z^(k) ,Z^'^ = p(Zj^ (k) |x(k)) ,p(£(k) IZ*^"^) / 

p(z^(k) IZ^"^) 

a 1 so f 

P(z2(*<)!Zj^('«)» 2.^''^^ • p(Z2(k)i'2^(k)) 

because of the indeoendence assumption. 

And sof 

D(x(k)|Z^) = (p(z„(k)!x(k)) / d(z (k)lz (k)fZ'^"^)]. 

- 2 - 2 1 

p(_x (k) I z ^(k) ,Z^~^ (a. 25) 

The terms on the right hand side of (^.25) are: 

p ( £( k ) 1 z^ ( k ) f Z^~^) --- has all the information about the 
states from processing up to the first measurement component 
z k ) , 

p ( Z 2 (k)|£(k)) --- is obtained directly from the 
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measurement equation and doesn't depend on the first meas- 



urement . 

o ( 2 ) I z ^ """ is obtained from the measurement 
equation but using all the results given by the processing 
of the first measurement component. 

This last term indicates that information is lost 
iff in a Extended Kalman Filter/ the linearization for pro- 
cessing a measurement is based on the predicted value 
instead of on the value estimated after processing the 
preceding measurement component. 

3 . The Tracking Problem 

In the tracking problem here addressed/ the measure- 
ments naturally occur at different times. Bearing indica- 
tions are basically continuous or, if preprocessed/ succes- 
sive measurements are available within short intervals of 
time. Frequency indications must be obtained from digital 
processing of records of considerable time length. Time 
delay indications are also obtained from digital processing 
but with longer time intervals. 

In earlier work/ 111/ [21 and [31/ these measure- 
ments are forced to occur simultaneously and are processed 
together . 



In the previous sections of this chapter it was 
shown that/ if the measurements are all available at the 
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same time/ they should be processed seoarately for better 
computing efficiency and/ with nonlinear measurements/ for 
possibly better estimation accuracy. 

Since one is now free from having to process the 
measurements together/ and encouraged to do so/ why not use 
the idle time of the computing eguipment and process the 
measurements at different times? 

If one does this the decrease in computation time 
obtained previously will be diminished because of the extra 
computation time reguired to do many predictions/ one for 
each time a measurement is to be processed. Nevertheless, 
the advantages are ootentially of great value, that is. 

--- a measurement can be processed as soon as it is 
available, with no waiting for other measurements. 

--- measurements that occur more often can be pro- 
cessed more times than measurements ocurring less fre- 
quent 1 y . 

--- the output of the filter is updated frequently. 

--- the one-step prediction is based on lineariza- 
tions through much smaller time-steps than before, hopefully 
with improved accuracy and fewer divergence problems. 

Figure 8 gives an idea of the improvement one 
expects to obtain by (1) processing simultaneous measure- 
ments separately and (2) processing the measurements as they 
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naturally occurs when a Extended Kalman Filter is beinq 
used . 



B. GRAPHICAL INTERPRETATION 

General understanding of the actions taken by a Kalman 
Filter during the processing of a measurement can be 
obtained through the graphical interpretation and visualiza- 
tion of these actions in a two-dimensional problem. 

From the conclusions of Section IV»A, it is advantageous 
to process any group of measurements with statistically 
independent noise components separately. In this section^ 
thus/ the processing of only one measurement at a time is 
cons i de red , 

1 , The Linear Measurement Case 

'tie will consioer a dynamic system with two state 
variables/ and a linear measurement of the form 

z(k) = H(k),)^(k) t v(k) = 

= hj^(k) (k) + h 2 (k),X 2 (k] + v(k) 

with v(k) a zero-mean/ discrete white Gaussian noise process 
with ECv^(k)l = r^(k). 

At time t a prediction is made from previous esti- 
mates and knowledge of the plant and system noise charac- 
teristics, Equations (4,3) and (4,4) are used for a linear 



83 





• I 







.processing simultaneous 
measurements separately. 



ih 






processing simultaneo 
measurements together 



us ^ 






-fi- 



\+l 




Figure 8 : Comparision of processing policies. 



systerDf (i.7) and (3.8) for a genera) nonlinear system. 

The Kalman equations (^.5)/ (^.6) and (^.7) are now 

applied to correct the predicted value according to the new 
measurement . 

Let's represent the quantities and relations 

involved as shown in Figures 9 and 10» where e and e are a 

~L “2 

pair of orthonormal basis vectors and a general vector is 
represented by 



^ = * 1^1 * * 2^2 



The following conventions and simplifications apply 



X* = X ( k ) 



X (k) 




X? 


1 




1 


X (k) 




*2 


-J 







X = X ( k ) = 




X 



1 



X* = x(klk-l) = 




~h, (kf 




"h " 


-1 


II 


1 


h, (k) 




h„ 


2 




2 
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Figure 9 : Graphical interpretation - Measurement 
lines and gradient. 



S6 




Figure 10 ; Graphical interpretation - Error ellipse. 
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I 




g = G ( k ) = 





II 




92 f ) 




^2 


- 




w — 



jh| = ^ ^2jL/2 



, , . , 2 ^ 2 , 1/2 
|gj - Ig^ + 92^ 



a= arc tan 

2 1 



e = arctan 



Let's define the function 



(J)(x_) = *^1*1 ^ h2X2 



(4.26) 



and call "measurement line" the locus of points where <j> is 
constant. Some of these lines are shown in Fig 9, 



Using this definition one sees that the previously 
defined h is the gradient of 6 with respect to x» or ^ f 
which in this linear measurement case is not a function of 
the states. 



The prediction error covariance matrix# 



P(k!k-1) = P' r 



11 



12 



12 



22 



is represented by an error ellipse whose direction ( 0 ) and 
axis dimensions (a - major axis# b - minor) are given by# 
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12] 



t an 


20 = 


= "12 ' 


Csl 
• CM 

a 

1 


) . 


. . -ir < 2 0 


a = 






p' / 
12 


s 1 n 


2 e 


b = 


‘'’ll 


♦ - 


p' / 
12 


s i n 


2 8 



a. The State Estimate 

Equation (^J.6) can be represented 
— ~ g»l'*esiduan 

where 

[residual] = z(k) - H ( k ) ,^( k I k- 1 ) = 

= (j) (£^*) v(k) 

and (4.30) means that a vector correction 
} g I . j res i dua 1 I is made to the predicted state 
tion of the vector q which has an angle B as 
ure 11. 



Equation (4.5) giveSf as shown in 



p ' h + p ' h 



11 1 



12 2 



p ’ h^ + 2p ' h h + p ' h^ + r 



11 1 



12 1 2 



22 2 



9 - 



"n ^1 



^ 12^1 ^ 22^2 
^ 2 Pl 2 ^^2 ^ ^ 22^2 ^ 



I (4.27) 

(4.28) 

(4.29) 

in the form 
(4.30) 

(4.31) 

of magnitude 
in the direc- 
shown in Fig- 

Append i X A ^ 
(4.32) 
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Figure 11 : Graphical interpretation - Gain vector. 



(1) The Direction of Correction 



The direction of g is seen to be indepen- 
dent of the measurement noise/ or/ given h and P(k!k-1) the 
direction of correction is already determined by 

e = arctan ^22^2^ ^ ^ ^12 ^2^^ 

(a. 33) 



From Equations (^.27) - (^.29) one obtains 



p' 

12 


- 1 
2 


( a - 


b ) . s i n 


26 




(9.3^0 


p ' 


- 1 


C(a 


+ b) + 


(a - b) .cos 


20] 


(a. 35) 


11 


2 








p ' 


- 1 


[ ( a 


+ b) - 


(a - b).cos 


20) 


(a. 36) 


22 


2 











Apolying these equations and the definition 
of a to t^,33) gi ves , 

(a/b - l).sin 2 e, + l(a/b + 1 ) - (a/b - 1 ) .cos 2 e) .tan a 

tan 3 = — ^ ; 

(a/b + 1 ) + (a/b - l).cos 20 + (a/b - 1 ) » s i n 2 9. t an a 

(a. 37) 

which can be reduced to 

( a/b ) . t an0 , ( t an 0 . t ana + 1 ) + (tana ” tan 0 j 

tan 3 = 

( a/b ) . [ t an 0 , t an a + 11 - tan0 .(tana ~ tan 0 ) 



(a. 38) 



91 



This equation shows that the direction of 



the correction generated by the Kalman Filter/ in this two- 
state oroblem, is a function only of the direction ( o ) of 
the gradient (^) and of the alignment (0) and the ratio 
(a/b) of the prediction error ellipse. 

Three special cases of interest occur when 
(i) a/b = 1/ (ii) b = 0/ and (iii)e=a. 

Case ( i ) - If a/b = t the error "ellipse" 
is actually a circle meaning that the prediction has the 
same uncertainty in any direction. No preferable direction 
pre-exists and the correction is in the direction of the 
gradient of the measurement function/ as can be seen from 
(a.i7) 

tan 8=10^ 2. tana) / (2 + 0) = tana 



or 



8 = a n T 



Case ( i i ) - If b = 0 the error ellipse is a 
line meaning that the prediction is exact along the minor 
axis and that only corrections along the major axis are 
necessary. As can be anticipated/ the direction of correc- 
tion will be the direction of the major axis/ as can be seen 
by taking the limit as b 0 of Eg. (^.30)/ which yields 

tan 8 - t an 6 
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m 





or 



3 — 0 +• n TT 

For any other intermediate value of a/b between 1 and » , 

the value of 3 will be between a and 9 ( t n n) , 

Case ( i i i ) - When 0 = a the gradient is 
colinear with the major axis of the prediction error ellipse 
and » f rom ( ^ . 38 ) » 

tan3 = ( ( a/b ) . t an 0. ( t an ^9 + 1) + 01 / ( ( a/b ) . ( t an^0 + 1 ) - OJ 
or 

tan3 = tan9 and 3 = 9+nir=a + n it 

Thus the direction of the correction is 
along the major axis of the ellipse which is colinear with 
the gradient of the measurement function. 

(2) The Amount of Correction 

The magnitude of the correction 

( I res i dua 1 I . I g ij ] can be better seen if one decomposes g into 
components m and n» along the gradient and the tangent to 
the measurement function^ as shown in Fig. 11. 

m = g^.cosa + g .sina (4.39) 

and from the definition of at 
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9 



sina = h 2 /|h_| 



(y .iJO) 



cos a = h 






Using Equations (^,32)/ (y.39) and (y.yO)r 



m 



p’ 

11 1 



+ p' h h 
12 1 2 



+ p’ h h 
12 12 



+ p ' h 2 
22 2 



p'h^ t2p'hh + p*h'^ 
11 1 12 1 2 22 2 



1 

lEf 



1 1 

= X . (y.yi) 

1 |hi 

where 

^ / (p* h^ + 2p' h h t p' h^) (y.y2) 

11 1 12 1 2 22 2 

Consider some special cases. How would the 
new estimate be computed if a noise-free measurement were 
obtained and the Kalman filter were aware of this fact» 
i . e . f r = 0 ? 



Note from Equations (4.^1) and (^.^2) that 
when the measurement is noise-freer r = 0, y = 0 and 

m = 1 / |Ji|| » a constant. 

No matter what P(kjk-l) is» the projection 
of the correction into the gradient is a constant if r = 0/ 
or t the tip of the correction vector follows a line normal 
to the gradient when P(k|k-l) is varied. 

Figure 12 shows the correction made to the 
predicted state when b = 0. In Appendix C it is shown that 
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for the present case of a linear measurement function the 



vector 0A» in the direction of h and length ( res i dua 1 1 / | h ^ 
has its end point A on the line 

<j>(x_) = 4>(x_‘) + residual - + (p (j^* ) + v - 



or 

'{'(x.) = V = 2 

which/ for r = 0 (i.e. v = 0) passes through the point x 

X.* . 



Since the correction has the direction of 
the prediction error ellipse (line) and has OA as its pro- 
jection into the vector h_» it can be seen from Figure 12 
that in this case of b = 0 and r = 0 the estimate will coin- 
cide with the true state. 

Figure 13 shows the correction for a gen- 
eral error ellipse but still with r = 0, It is seen that the 
projection into h is the same as before. The direction of 
correction 0 is somewhere between a and 6 . 

A general correction for a noisy measure- 
ment is shown in Figure 1^/ with different scaling than the 
previous figures. The direction of correction g is first 
determined and suppose it is as shown. The projection of the 
correction into h^ would be OA if r were 2 eroJ for r 0 the 
reduction factor is applied and suppose the projection is 
vector OB as shown in the figure. Taking a normal to from 
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Figure 12 : Linear measurement - Perfect estimation 
for b=0 and r=0 . 
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Figure 13 : Linear measurement - r=0. 
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point the point C where it cuts the direction of correc- 
tion determines the final correction OC and the new esti- 
mate. 



b. The Estimation Error Covariance Matrix 

The effects on the estimation error covariance 
matrix can be seen in the following way. From Equation 
(a.7) , 

P(k) = II - G(k)H(k)lP(k|k-l) 



or 





P-. o 




1 - g h. 


• d h “] 






p ' 


11 


12 




1 1 


1 2 




11 


12 


P 


P 




- g h 


1 - g h 




P • 


P’ 


12 


22_^ 




2 1 


2 2 




12 


22 



















Using the values of g^^ and P 2 from (4.i2) one 



finds that 



2 2 2 
Ch^I P j_ 3_ ^22 * ^12 ^ ^ ^ 7 c 



^ 12 = ■ ^ 11 ^ 22 ^ " ^12 " 



p = th^(p* p* - p’^ ) + p’ r2) / c 

22 1 11 22 12 22 



(a. 43a) 
(4. 43b) 
(4.43c) 



where 



c = I p,*, h^ + ?D* hh + p* h^'*-r^l 
11 1 12 1 2 22 2 



(4.43d) 



98 




0 ^ 

/ 

f 

/ 

t 

I 

Figure 14 ; Linear measurement - General correction. 
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I 




And the estimation error ellipse has the orien- 



tation 



tan 26 = 12 ^ 



- P22 ^ 



(b! 



2hih2 (p’j_2 



- P 



11 "^22 



2p ' 

^12 



- - 



"u "22 > » 



(P' 

11 



• Poo 
22 



(a.aa) 



The limiting situations occur in the special 
cases when (i) r -> 00 and (ii) r = 0, 

Case ( i ) - When r is relatively larger little 
information is provided to the filter; the estimate is 
essentially the same as the prediction^ and the estimation 
error ellipse is essentially the same as the prediction 
error ellipse. This can be seen from Equations (y.^3a) 
(y.y3d)/ as r-voo. 

Case ( i i ) - when r = Of all errors in the pred- 
iction/ in the direction of the gradient of the measurement 
function/ are corrected and the estimation error ellipse 
becomes a line colinear with the measurement line. As shown 
below/ this result is obtained from Equations (^.^3) and 
when r is made equal zero. 

From (4.a3a) - (A.^3d)/ p .0 is easily shown 

11 22 

2 

to be equal to p^^ r = 0/ indicating that P(k) becomes 

singular. From (^1.^4)/ 

tan 26 = (h^ ' ^ 2 ^ = 
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2a 



= 2(h /h ) / tl - (h /h )2 ] = 

2 1 2 1 

2 

= 2 tan a / (1 - tan a) = tan 

From analysis of the individual terms of (^.^ 3 a) 
- (^I.U 3 d)» this result is to be interpreted as 

2 6 = 2 a +. 

or 

6 = a ^ ■nr/2 

When the minor axis of the 
already zero^ as in Figure 12^ i.e.f 
measurement is noiseless^ then P(k) = 0 and complete cer- 
tainty about the system state is obtained. 

These two last paragraphs show a very interest- 
ing situation: if two perfect measurements are made on a 
two~state plant and the filter is aware of this (i.e. r = 0 
is used in the Kalman equations)# then an exact estimate is 
obtained --- the estimation error covariance matrix col- 
lapses into the zero matrix. 

For a general situation with normal measurement 
noise values# the estimation error ellipse is rotated toward 
the measurement line# away from the gradient# with 6 occu- 
pying a position between 9 and a + tt/2. 



prediction ellipse is 
^12 " ’^1 ^^22 ' 
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The general expressions for the lengths of the 
major and the minor axis of the estimation error ellipse can 
be obtained by the application of Equations (^.28) and 
(a. 29) to (a.asa) - (9.a3c). 

c. The Important Points 

Summarizing the principal points seen up to now» 

--- the direction of correction in the estimate is 
defined at the end of the prediction phase and doesn't 
depend on the measurements. 

--- if the measurement is noiseless and r = 0 is used 
in- the Kalman filter equations# the projection of the 
correction on the gradient of the measurement function is 
(residual / fJlp' independent of PCklK-l), The estimation 
error covariance matrix becomes singular indicating# in a 
two-dimensional problem# that the estimation error ellipse 
becomes a line --- a point when the prediction error ellipse 
is al ready a line. 

--- if the measurement is noisy# the projection of the 

correction on the gradient is multiplied by the reducing 

2 

factor ll / (1 + Y error ellipse is rotated toward 

the measurement line# away from the gradient of the measure- 
ment function. 

--- if the measurement is too noisy# the reducing fac- 
tor tends to zero and the estimation is the same as the 
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I 




predi ct i on 



2 . The Nonlinear Measurement Case 

Assume the same two-dimensional system as before but 
with nonlinear observations of the form 

2 (k) = h(x(k),t ) + v(k) 

— K 

a* The Extended Kalman Filter Approach 

The prediction phase is the same as before but 
the estimation equations are now (3,18)» (3.19) and (3.21) 

for an Extended Kalman Filter. 

Let's, again define a measurement function by 

, 4>(^) = h(^) 

and redefine the gradient of 6 with respect to £. 

Now the components of ji are no longer constants 
but depend on ^L* the predicted point these components 

will be 

h = A h(;(k!k-n ) , h hC»(k;k-l)) 

3*1 - 3X2 “ 

Equations (^.27) through (^.4A) are all still 

va 1 id. 



Let's consider two tyoical measurement functions 
and observe graohically how the new estimate is obtained. 
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Case ( i ) 



assume the measurement function is of 



the form 

(}i(j<) = a ret an C(x^ ” ^ ^ ^*2 ” ^2 ^ ^ (^.^5) 

The measurement lines will be as shown in Figure 15. 

Suppose now a noiseless measurement is obtained 
and is to be processed. Fig. 16 shows the situation when the 
prediction error ellipse is a line (b = 0). In Appendix C it 
is shown that. for this special measurement function the vec- 
tor OA of length t res i dua 11 / | h | will end somewhere before a 
point B on the curve (j)(£) = z. Supposing its length is as 
shown/ it can be seen that an imperfect estimate is 
obtained. Nevertheless the estimation error covariance 
matrix/ given also by Equations (4.^3)/ will collapse into 
the zero matrix indicating/ falsely/ that a perfect estimate 
was obtained. 

Figure 17 shows a typical situation for a gen- 
eral prediction error and a noise-free measurement. The tip 
of the correction vector can be seen to follow the line A - 
A' parallel to the line (|)(j^) = ((>(£') instead of being along 
the line = ({)(£*). As seen in the previous section/ 

the estimation error ellipse will erroneously be determined 
as a line also along A - A', indicating falsely the direc- 
tion of the true state. 

Case (ii) - assume the measurement function is 



of the form 




Figure 15 ; Nonlinear measurement lines I. 



105 




i 



Figure 16 : Nonlinear measurement I - Imperfect 
correction for b=0 and r=0 . 




Figure 17 : Nonlinear measurement I - General 
correction for r=0. 
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(a.ab) 



The measurement lines are as shown in Figure 18. This same 
figure shows the new estimate obtained when the measurement 
is noiseless (r = 0) and the prediction error ellipse is a 
line (b = 0), The same observations of Case (i) are valid. 
As shown in Appendix C the point A is on the curve ,j, ( x ) = z 
for this special measurement function. 

These figures have their values mostly because 
they show» in a simple way» the large errors that can be 
made in the processing of nonlinear measurements by an 
Extended Kalman Filter. Two types of errors are generated; 
one in the determination of the new state estimate# the 
other is seen in the fact that the estimation error covari- 
ance matrix P becomes a poor approximation to the true error 
covariance# making the filter non-optimal. These difficul- 
ties are a function of the nonlinearities involved and# 
apparently# are more significant when the predicted values 
are poor and when the measurement noise is assumed to be 
small when evaluating the extended Kalman filter equations. 

b. The Iterative Approach 

A way of correcting# in part# these deficiencies 
is the adoption of iterative procedures. The first improve- 
ments can be obtained if# after the Extended Kalman Filter 
has been applied to correct the prediction# the filtering 
process is repeated but with the linearizations made about 



- a, ) t ( 



- a_ ) 
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Figure 18 



Nonlinear measurement II - Imperfect 
correction for b=0 and r=0. 



the just obtained estimate# hopefully a better approximation 



to the true state than the crude prediction. 

Analytically the process consists of replacing 
the estimation equations (3,18)# (3.19) and (3,21) by the 

iterative equations 

/ 

i = _x(k|k-l) + 6 l 2 (k) - _h(i(k|k-l))] (4,47) 

i+1 i 



= P(k;k-l)Hj^(K P(k!k-l)h^ + 


R) ^ (4.48) 



H . = -2- h(x.) # i=0# 

■ ^ 3x 1 

where 


•••# f 


X =x(k!k-l) # 

— o — 


1 X > 

II 

1 X> 
<— > 
X" 


P(k) = II - G^H^)P(klk-l) 


(4.49) 


For our two-dimensional 
can be expressed as 


problem Equation (4,47) 


X = X* + q .(residual! 

-i+1 - 

where 


O 

in 

• 


residual = z - h(x') = <f> (_x*) • 


• (x • ) + v(k) (4.51) 



Adding and subtracting cji ( )^ ) to the residual# 

i 
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or 



X = X * 

~i+l - 






<()(x ) + (j> (x ) 

^ ~1 —1 



(|) (x • )1 



x\_^ = X • + Iz - d) (x . )1 - q. . l(|)(x ’ ) - (|) (i . )] (4,52) 

These equations can be applied until there is no 
significant difference between consecutive estimates/ and 2 
or i iterations are normally enough. 

The effect of these equations can be easily 

visualized for the simple case of b = 0 and r = 0, The 

sequence of Figures 19 - 21 show this effect for one of the 

special measurement functions already studied. In Figure 19 

the first estimate is obtained as done previously in Figure 

16, The vector OA/ which is smaller than OB as shown in 

Appendix C/ represents the orojection of the correction on 

the gradient h , A normal to h from point A cuts the pred- 

— o 

iction error line (b = 0) at the first estimate x^. 

In Figure 20 the terms of Equation (4,52) are 
shown for i = 1, Two correction terms are applied to the 
prediction )^'/ both in the direction of g_^ but with dif- 
ferent residuals. The direction of g^ is again along the 
prediction error line (b = 0), 

The gradient of the measurement function at the 
first estimate The projection of the first correc- 

tion ( 2 .j^(z - d>(x.^)l) into h^^ is smaller than CE/ say CO, 
Taking the normal at D one gets the first correction term 

1 1 1 




Figure 19 



Iterative process - First estimate. 
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CH 



The projection of the second correction term, 

( t ))), into h is smaller than CG, say CF, 

i - i —1 

The normal from point F determines the second correction 
term ( IC) . 



Adding the three vectors CH and IC one has 
the second estimate, shown in Figure 21, and the next itera- 
tion can be processed. 



If in Equation (^1.^7) one modifies the residual 
by 

Z_( k ) - ]^(^(k|k-l)) = ^(k) - fil(x. ) + 

i 

+ H ,(x(k|k-l) - X ) + ...] = 

i — — i 

= z(k) - h(J . ) - H . [i(kjk-l) - X ] 

— — —1 1 — — i 

then Equation (U,il7) becomes 

i , = x(klk-l) + G l 2 (k) - h(x ) - H tx(k|k-l) - x ]] 

“1+1 — i — i i“ “1 

(a. 53) 



This variation in the residual is very small to 
affect considerably the estimates but in 111) it is shown 
that its use produces a better agreement between the calcu- 
lated estimation error covariance matrix and the true error 
covariance, when processing a noisy measurement. 
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Figure 20 : Iterative process - Correction terms. 
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Figure 21 : Iterative process - Second estimate. 
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3. Summary 



The direct application of these graphical interpre- 
tations to a general problem is impractical. Nevertheless 
the basic points observed here can be generalized for Kalman 
Filters: 

--- the "direction" of correction in the estimation 
is defined by the predicted quantities and doesn't depend on 
the measurements. 

--- the amount of correction in an estimation^ along 
the direction defined^ is the result of a weighting process 
between the uncertainty in the measurement and the confi- 
dence in the predicted values. This correction is a function 

2 

of the factor [1 / (1 +-Y )J where/ in general/ for a n- 

dimensional system 

n n 

o’., h.h, (a. 51) 

x-1 j=l xj 1 j 

--- the errors in the estimate from the processing 
of a nonlinear measurement can be large and/ most imoortant/ 
the estimation error covariance matrix may become a' very 
poor approximation of the true error covariance/ forcing the 
filter to inaccurately weight future measurements. This 
problem depends on the nonlinearities involved and seems to 
be worse when the predicted values are poor and/or the meas- 
urements are considered by the filter to have relatively low 
noise. 
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the errors indicated above may be reduced with 



the use of the iterative equations (^.48), (4.U9) and (4.47) 
or (4.53)> instead of the standard approach given by the 
Equations (3.l8)f(3.19) and (3.21)^ at the expense of 
greater computing effort. 



117 



if 










¥ 



av 




V. COMPUTER SIMULATION 



A. GENERAL CONSIDERATIONS 

A conoputer simulation was implemented in a Unix Operat- 
ing System running on a POP 11/50 with the FP IIB floating 
point processor^ located at room 506 of Spanagel Hall/ NPS. 
The whole program is interactive and is expected to be 
self-explanatory to a user. Its block diagrams and coding 
are presented in Appendix E, 

The graphical outputs were generated on a Tektronix 
yOl^J-l terminal. The presentation of the results of the 
various simulations is simplified by the adoption of the 
following conventions. 

The metric system was used in all calculations and out- 
puts. The axes are dimensioned in kilometers. Heading and 
bearing angles/ expressed in radians/ are measured counter- 
clockwise from the positive X-axis. 

' Figure 22 shows a typical plot. Points of the true track 
are represented by the sysmbol and interconnected by a 
solid line. Buoys are represented by a distinctive boxed 
symbol. The portion of the true track that is in the range 
of a buoy is bounded/ if necessary/ by dash-dotted lines. 
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Figure 22 : _ Tru e an d e s t imated traj ect ories 
A typical run. 
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The signals received by the sonobuoys are assumed to be 
sent to a ship or aircraft where they are processed by spe- 
cialized equipment to generate certain measurements. The 
major characteristics of these measurements are printed on 
top of each figure and obey a code whose general form is 

L N1 N2 N3 S 

where 

L - a capital letter indicating a type of measurement. 
B stands for bearing? F for frequency? D for frequency 
difference? R for frequency ratio? and T for time delay. 

N1 - an optional number indicating^ when necessary^ the 
buoy that generates the measurement. If this number has two 
digits it identifies two buoys^ thuSf T 12 indicates a time 
delay measurement between buoys 1 and 2. 

N2 - a number indicating how many seconds after the 
start of the problem this measurement will be first 
obtained. It is also optional. 

N3 - a number indicating the period between consecutive 
measurements. 

S - a number indicating the standard deviation of the 
simulated zero-mean Gaussian noise that is associated with 
this measurement/ followed by an optional small letter indi- 
cating dimension: d for degrees, h for Hertz and s for 
seconds . 
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The interpretation of the code presented in Figure 22 is 
then that both buoys 1 and 2 provide bearing and frequency 
measurements with a period of 100 seconds between consecu- 
tive set of measurements. Buoy 1 provides its first set of 
measurements 100 seconds after the start of the problem? 
buoy 2» perhaps because of its limited ranges provides its 
first set of measurements only after 500 seconds of the 
start of the oroblem. The standard deviation of the measure- 
ment noise is 5 degrees for the bearings and 0,0^ Hertz for 
the frequencies. 

The filter receives the measurements for processing. It 
is given some initial conditions and the initial position is 
marked in the figures by the first diamond with the letters 
I . P . on top. 

As the target moves and the measurements are simulated^ 
the filter generates estimates of the target parameters 
which are represented by small diamonds/ interconnected by 
dashed lines to represent the estimated path. The diamonds 
and the symbols correspond to the same time points. 

Associated with each estimated point the filter provides 
an estimation error covariance matrix which indicates the 
amount of confidence or/ conversely/ uncertainty that should 
be given to this estimate. An interesting way to show this 
in a trajectory plot like Figure 22 is by plotting the one- 
sigma error ellipses associated with each estimated point. 
These error ellipses have their orientation and axis dimen- 
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sions given by Equations (^,27) - (^,29), and can be seen in 
Figure 23 . 

These figures were obtained by running the simulation 
once. They may represent a typical real-time result. 
Nevertheless they do not give the expected or average 
behavior of the filter to the situation simulated. To obtain 
this average behavior many Monte Carlo runs have to be exe- 
cuted and their individual results processed to generate the 

/ 

sample statistics. The sample statistics formulas for 
obtaining the sample means and covariances can be founds for 
example^ in page 86 of Ref. [21. 

In this work all the statistical results were obtained 
from 200 MonteCarlo runs. Figure 29 shows the average 
behavior of the filter after these runs. Now the small dia- 
monds represent the sample mean estimated position and the 
ellipses show the spreading of the results about the mean 
values. Note that at the initial position oointr since in 
all runs it was the same» the ellipse is the point itself. 

To observe the errors in the estimation of the speed and 
heading of the target different plots have to be used. Fig- 
ure 25 shows the sample mean error in the estimation of 
speed as a function of time. Figure 26 shows the sample 
mean error in the estimation of the heading of the target. 

Other types of plots or outputs could be generated from 
the results given by the program simulation» like the sample 
covariance between estimated velocity and bearing» if it 
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ever becomes necessary 



B. SELECTED RESULTS 

As has been seen in this work» the number of variables 
in the passive tracking problem is very large. Some specify 
the true target track and the true position and characteris- 
tics of the sonobuoys? others define the measurement 
scheduling and noise characteristics. The nonlinear filter 
that processes the measurements has its own variables and 
must also assume values for the other variables^ which may 
be different from the true ones. 

For any combination of these variables a new problem is 
formed and its simulation may generate a number of interest- 
ing conclusions. The results presented are reoresentative of 
those obtained from the execution of the simulation program 
described in Appendix E. In most of the results presented a 
certain reasonable initial condition for the filter was 
assumed? the position of the buoys was assumed constant and 
known? the ranges within which reliable bearing^ frequency 
and time delay measurements could be obtained were assumed 
edua 1 . 

1 . Non-maneuvering Target 

a. One Buoy in Action at a Time 

Figure 27 shows the behavior of the filter in a 
situation where only one buoy is in contact with the target 
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at any instant. The range and position of the buoys were 
pre-arranged to satisfy this reauirement. Target soeed is 

5.0 m/s; its heading is 315 degrees measured from the posi- 
tive X-axis. The initial assumptions of the filter include 

6.0 m/s for speed and 325 degrees for heading. The measure- 
ments were processed as suggested in 121. Frequency and 
bearing measurements are available every 100 seconds and 
processed simultaneously; standard deviations of 5 degrees 
and 0.0^4 Hert? were assumed for the noise components. 

The frequency measurement provides^ indirectly^ 
some range information and» this way/ complements the bear- 
ing measurement. This range information can also be obtained 
by judiciously positioning the next buoys or by having two 
or more buoys in contact with the target at one time. 

In Figure 20 the frequency measurement was elim- 
inated and the simple bearing measurement allowed the esti- 
mator to track the non-maneuvering target almost as well as 
before/ because of the special position of the buoys. Note 
that the ellipses tend to align with the measurement lines. 
This can be explained from the discussions of Section IV/B 
and Appendi x D . 

There are many ways to improve this tracking. A 
very small improvement can be obtained if less noisy or more 
frequent frequency measurements could be obtained. To get 
more precise frequency indication one needs better resolu- 
tion from the FFT processors; to obtain better resolution 
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Figure__ 27 ; .True and estimated .tralectories f or a 
non-maneuvering target - Simultaneous 
frequency and bearing measurements. 
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one needs longer records which/ 



in turn/ introduce more 



errors due to target -sensor geometry variations during 
record time. To reduce the period between consecutive 
independent frequency indications one has to sacrifice reso- 
lution but this would reduce even more the effect of the 
frequency measurement on tracking accuracy. Besides/ an 
emitted frequency of 300 Hertz from a target with a relative 
velocity of 5 knots toward a buoy is shifted by only about 
0,5 Hertz and it can be seen that not much resolution can be 
spared. 



The bearing information is basically continuous 
and it seems reasonable to admit that only economical limi- 
tations exist for obtaining less noisy or/ equivalently/ 
more frequent bearing measurements. Figure 29 shows the 
effect of processing only bearing measurements/ but more 
frequently (at 25-second intervals). In Figure 30 the stan- 
dard deviation of the measurement noise was reduced to 1 
degree. Note that the alignment of the ellipses with the 
measurement lines is more pronounced with less noisy meas- 
urements. Mote also that the error in range while the first 
buoy is in contact with the target is not corrected by 
improving the accuracy of the bearing measurements. 

Figures 31 and 32 show the mean errors in 
estimating the velocity and bearing of the target. The sym- 
bol is associated with Fig 27/ the triangles with Fig 28 
and the diamonds with Fig 29, 
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Figure ^ 29 __ •. .True_ and estimated, trajectories for ,a_. 

non-maneuvering target - Reduced inter 
val between consecutive measurements. 
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Since the contribution of the frequency measure- 
ments is so small compared to that of the bearing indica- 
tions in this case f the partitioning of the measurements has 
very little effect on tracking accuracy. Also the iterative 
techniques suggested in Section IV»B,2 have no marked effect 
on the processing of the frequency measurements for the same 
reason . 



b. Two Buoys in Action at a Time 

The placement of another buoy in action can 
greatly improve the tracking accuracy^ as shown in the next 
figures. The buoys are now arranged so that two buoys^ 
instead of one^ are in contact with the target at any time. 
The scaling is changed so that the estimation errors may be 
better appreciated. 

In Figure 33 the added buoy is a Lofar buoy that 
can only provide frequency measurements. The two frequency 
and one bearing measurements available every 100 seconds 
were processed simultaneously. Note the improvement in accu- 
racy when buoys 3 and A gain contact with the target in 
replacement of buoys 1 and 2. This also is explained from 
the discussions of Section IV^B. 

In Figure 3^1 all buoys contribute with frequency 
and bearing measurements which are still processed simul- 
taneously every 100 seconds. The position of the buoys^ in a 
line perpendicular to the true track, is not a very good 
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Zigure .33, ..:.,True ,and._estimated . trajectories f or a 
non-maneuvering target - Simultaneous 
processing of one bearing and two fre 
quency measurements. 
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position for two Difar buoys. Nevertheless the resulting 
tracking accuracy is good for this non-maneuvering target. 

In Figure 35 the two frequency measurements/ 
which require a five-state plant for their processing/ were 
replaced by a frequency difference measurement with 0.0^1 
Hertz of standard deviation in its noise component/ which 
requires only a four-state plant/ as discussed in Section 
11/0/2. The results are basically the same/ mainly because 
they are mostly dependent on the bearing measurement and not 
on the frequency information. The use of the frequency ratio 
measurement also gave the same kind of results when a very 
low noise measurement/ with about 0.0001 units of standard 
deviation/ was considered. 

In Figure 36 the frequency measurements were 
eliminated and only bearing measurements were processed. 
Because of the bad position of the buoys/ on a line perpen- 
dicular to the true track/ the accuracy of the tracking 
along that line is not good. The ellipses show a reasonable 
spreading of the estimates along perpendiculars to the true 
track. Figure 37 shows the same plot/ with different scal- 
ing/ for the first 1000 seconds of the path. 

For different positions of the buoyS/ especially 
if one moves the buoys off the perpendicular to the true 
track/ better tracking can be obtained as shown in Figure 
30. The tendency of the ellipses to align with the measure- 
ment lines explains the improvement obtained. 
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Figure 35 : True and estimjated_ trajectories for a 
non-maneuvering target - Simultaneous 
processing of bearing and frequency- 
difference measurements. 
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non-maneuvering target - Separate pro 
cessing of bearing measurements - 
Buoys in different position. 
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By reducing the period between consecutive bear- 
ing measurements and processing them separately one may 
obtain the accuracy shown in Figure 39. The ellipses are 
like the ones of Figure 38 but with reduced dimensions. 

Good results can also be obtained if one 
includes the time delay measurements. Figure 90 shows the 
result when they are applied to this situation. Notice that 
the ellipses are very small reflecting the fact that very 
good information about the position of the target is given 
by a low-noise time delay measurement. Also the ellipses 
tend to align with the tangent to the measurement lines 
defined in Section IV»B which^ in this case of time delay 
measurements; are hyperbolas. 

2 . Maneuvering Target 

The next paragraphs discuss results obtained when a 
selected target maneuver is simulated --- a total turn of 
180 degrees at 9 degrees per minute with a constant speed of 
5.0 m/s . 



a. One Buoy in Action at a Time 

Figure 91 shows the tracking obtained by pro- 
cessing frequency and bearing measurements provided by a 
single buoy located close to the center of curvature of the 
target path. The filter takes some time to react to the 
target maneuver and; when it does» an ove r-co r rec t i on is 
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applied. The result is that the estimated path is very dif« 
ferent from the true oath of the target. About the same 
result is obtained when the buoy is placed on the other side 
of the tracks i.e.» a small delay to react to the maneuver 
followed by an over-cor rect ion . 

In Figure ^2 the frequency measurements were 
eliminated and only bearing indications were processed to 
generate an estimated track very close to the one of Figure 
^l. The alignment of the ellipses with the measurement 
lines is again very clear and may suggest where two buoys 
should be placed in order to improve the tracking. 

b» Two Buoys in Action at a Time 

Figure ^3 shows what can be obtained by process- 
ing the measurements provided by two buoys. The range of the 
buoys was adjusted so that both are in contact with the tar- 
get during all the path shown. Both are Difar buoys and pro- 
vide bearing and frequency measurements with a period of 100 
seconds^ which are processed simultaneously by the filter. 

In Figure 0^1 the measurements were processed 
separately and a small improvement in the tracking was 
obtained at the end of the path/ although during the 
maneuver the simultaneous processing worked better. It was 
observed in other simulations that the degree of improvement 
obtained by partitioning the bearing and the time delay 
measurements is dependent on the maneuver of the target and 
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the relative position of the buoys. The worst situation is 
when the target maneuvers toward the buoys. This happens 
when the partial estimates are worse than the predicted 
values and occur sometimes during fast maneuversr or slow 
maneuvers with long intervals between consecutive measure- 
ments . 

In Figure ^5 the low-noise time delay measure- 
ment was added and» instead of reduc ing the errorSf the 
tracking became worse. The direction of the ellipses is 
clearly along the tangent to the hyperbolas of position^ 
indicating that the filter is uncertain about the range of 

the target to the buoys during all the maneuver. This 

/ 

characteristic of the measurement lines in this case/ and 
the fact that the target is modelled by the filter as fol- 
lowing basically a straight path/ explain the form of the 
estimated track. As soon as the target reassumes a constant 
path the filter starts to recover. 

In Figures 46 and 47 the iterative techniques 
described in Section IV/B/2/b were applied to the frequency 
and time delay measurements. In Figure 46 only one itera- 
tion was executed; in Figure 47 the iterated formulas were 
applied three times. The improvement is clearly seen by com- 
paring with Figure 45. 

As with the partitioning of the measurements/ it 
was observed from other similar simulations that the appli- 
cation of the iterative techniques does not always provide 
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better tracking. It depends on the type and direction of the 
maneuvers/ the position of the buoys and/ consequently/ on 
the position of the measurement lines. 

In the case of the non-maneuvering target the 
use of two buoys instead of only one was almost always suf- 
ficient to provide good tracking. If for the maneuvering 
target one uses three buoys judiciously positioned/ instead 
of only one or two/ one can also obtain good tracking as 
clearly shown in Figure ^<8. 



C. SUMMARY 

Below are listed the most important facts observed from 
the simulations: 

--- The tendency of the error ellipses to align with the 
measurement lines was observed in Section IV/0. The simula- 
tions show that the sample covariance ellipses also follow 
this trend. 

--- The tendency of the error ellipses to align with the 
measurement lines can be of great help in the practical 
solution of filtering situations. If/ for example/ the error 
ellipse is very thin/ i.e./ have a very high ratio a/b/ then 
the best measurement to process is the one whose measurement 
line is approximately perpendicular to the major axis of the 
ellipse. In the tracking oroblem this can be obtained by a 
bearing measurement from a buoy placed along the perpendicu- 
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lar to the major axis of the ellipse^ or by a time delay 
measurement from two buoys placed along the direction of the 
major axis of the ellipse. 

--- The frequency measurement have a very small effect 
on thetracking accuracy. 

--- The use of frequency difference or frequency ratio 
measurements improve computing efficiency (by eliminating 
the need for having the rest frequency as a state in the 
filter model) without noticeable effect on accuracy. 

--- The partitioning of the measurements improves com- 
puting efficiency/ may provide better tracking accuracy and 
allows the processing of measurements as they occur/ and 
thus is of great practical importance. 

--- For low noise measurements the iterative techniques/ 
although requiring an increase in computing/ are capable of 
providing better tracking for maneuvering targets. 

--- Tracking with only one buoy is acceptable only for 
non-maneuvering targets. 

--- The relative position of the buoys is one of the 
most important factors which influences the quality of the 
t rac king. 
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VI. CONCLUSIONS 



A, SUMMARY OF RESULTS 

The optimal estimation of characteristics and/or parame- 
ters of a certain class of nonlinear/ dynamic/ stochastic 
systems has been studied in a probabilistic environment 
using Bayes formulation concepts. Approximate solutions and 
filtering algorithms were generated and/ among them/ the 
known Extended Kalman Filter equations and higher order 
filtering equations have been obtained through this method. 

The problem of tracking submarine targets using special 
passive sonobuoys was modelled and with this model extensive 
simulations were executed to allow the study of the problem 
in detail. 

Most of the results indicate that the frequency measure- 
ments have minimal effect on the filtering process. The 
small contribution to range information that thev do pro- 
vide/ when associated with bearing measurements/ can nor- 
mally be obtained otherwise by Judicious placement of subse- 
quent buoys. 

The utilization of other types of measurements Such as 
the frequency difference/ the frequency ratio and the time 
delay measurements/ proves to be a great help in improving 
computing efficiency by eliminating the rest frequency as a 
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necessary states and thus reducing the dimensionality by 
one. Also» especially with the time delay measurements a 
great improvement in tracking accuracy is possible. 

The concept of partitioning the measurements and pro- 
cessing them separatelys even if they occur s i mu 1 t aneous 1 y r 
is shown to bring advantages in computing efficiency and 
alsos for nonlinear measurements/ in tracking accuracy. 
This concept is of great practical imoortance/ for» in 
situations where the measurements are sporadically and non- 
periodically received it is most important to be able to 
process the measurements as they naturally occur. 

The graphical interpretation of the action of Kalman 
filters/ developed in this work/ provides insight into the 
importance of each variable of the problem in the . filtering 
process. The direction and magnitude of the correction which 
is applied to the predicted values to generate the new esti- 
mate values can now be anticipated as a function of the 
measurement . 

Nonlinearity errors have been graphically presented and 
iterative techniques/ including the known Iterated Extended 
Kalman Filter equations/ have been suggested to counteract 
their disruptive effect. The application of these tech- 
niques to the tracking problem shows that improvement in 
tracking accuracy is possible. 

The graohical interpretation also indicates the very 
practical conclusion that the error ellipses tend to align 

I 
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with the measurement linesf as defined in Chapter IV. This 
provides guidance for optimal positioning of the buoys and 
the types of measurements to process. This observation is 
reinforced by the simulation results of Chapter V, 



B, SUGGESTIONS FOR FUTURE RESEARCH 

The concept of partitioning the measurements allows one 
to. process the measurements separately and opens some ques- 
tions for future research^ such aS/ in which order should 
nonlinear measurements be processed to obtain maximum effi- 
c i ency? 

The graphical interpretation of the filtering process 
allows anticipation of the direction and magnitude of the 
corrections which are then applied to predictions to gen- 
erate new estimates. This situation suggests future research 
in the determination of the optimum characteristics of the 
measurement functions which in turn can determine the 
optimum positions and characteristics of sensors. 

The expansion of the model to include accelerations 
should be considered if maneuvering targets have to be 
tracked efficiently and the increase in computing power can 
be obtained. Consideration of the depth of the target and 
the inclusion of uncertainty about the position of the buoys 
are additional ways to extend this study. 
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The model and simulation implemented herein/ 



w i t hout 



benefit of classified information/ is very flexible/ but 
many idealizing assumptions have been made. IhuS/ a natural 
extension of this work is to apoly the algorithms and con- 
cepts developed to a real problem using actual sensor data. 
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APPENDIX A : PROBABILITY RELATIONS 



The relation 



X 



a + G . V 



(A. 1) 



is given where x and v are random vectors of dimensions n 
and pf respect ivelyr ± is a n-dimensional deterministic vec- 
tor and G is a (nxq) matrix of deterministic coefficients. 

The joint density is given and the joint density 

p^(x.) is wanted. 

Case ( i ) - If n = q and G ^ exists/ the solution is 
given by [131/ 



V = G'^x - al = f(x) 



(A. 2) 



and 




(A. 3) 



Case ( i i ) - If n < q one could add dummy variables 




and create 



x* = a’ ♦ G' .V 
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3 



where 



X 



I 



X 



rril 



X 

<U 



a' 






Now f if G'~l exists^ 

V = lx' - a’l = f(x») 



From Case (i)/ the solution is 



P ,lx’ ) = 

X 







. P (f (X ' ) ) 
V — — 



and 



p ( X ) = 

X • ~ 



p ,(x ' ) .dx , ,.dx . ,dx 
X “ n+1 n+2 q 



Case ( i i i ) - If n > q the following 
app lies. 

Let 



(A. a) 



(A. 5) 



(A. 6) 



de ve 1 opmen t 



16 a 



r X “I 
*1 
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*2 




L.' 
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• 






9 _a — 


... 


q +1 

• 




q+1 

• 




\+i 

• 


X 




X 




a 


n 




L " - 




n 



and 



From 

X.' = 

t hen f 

V = 



G* 




G* 


qFl,l q+i,q 




a 

^+1 

• • • 


^ 

n , 1 ri 5 q 






L ^ 




L J 



(A.l), one also has that 
a' ♦ G' .V 
if G' ^ existSf 
G C X ' - a ' ] = f ( X • ) 



(A. 7) 



From Case (i) one then has 






. f X ) 

q 







. p <f («' ) ) 

V — — 



(A. 8) 



but 



Now consider the variable x ... Given (A,7)r 

q+1 



p(x,, !x.X .,..fX ) = o(x !v»v»...fV ) 

q+1 12 q q+1 12 q 



X ^ Q V V +♦••+0 V 

q+1 q+1 Vl»l-1 q«-,2 2 ^q+l,q q 



= a + g .V 
q+1 ^+1 - 



or 
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-q+1 ■ ®q+l %+l *-- 



. f (x ' ) = 1 



q+1 



(A. 9) 



From (A. 9) comes 



q+1 12 q 



q+1 q+1 



and 



p(x fX /.../X ) = 6(x - 1 ),p (x tx ,.,.,x ) 

1 2 q+1 q+1 q+1 x* 1 2 q 



Doing the same thing for x ...f x , one finally 

q+2 n 



gets 



n 



p (x) 

X •— 



6 ( X 



1 . ) 
1 




detl -9_f 
axT 



p (f(x’)) 

V — — 



i=q+l 



(A. 10) 



For the special case of q = 1/ the vectors _x ' / a’ and 
g' become scalars and the sbove equations lead to 



f(x') = f(x) = (x ~ a ) / g 
1 1 11 . 






1 = a + (x - a )g /g 

r r 1 1 r 1 



and Equation (A. 10) becomes 
n , 



p (x) = 

X — 



6(x. - 1 . ) . 

11 I 1 



. p ((x -a )/g ) 
V 1 11 



i=2 



(A. 1 1 ) 



Another similar formula can be obtained for this soe- 
cial q = 1 situation. The matrix G in (A.l) is now a n- 
dimensional vector and/ if one breaks (A.l) into its n 
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individual scalar equations one can obtain 



p(x^!v) = 6(x^ ’ 



Also note that 

p(x^|Xj,v)=p(x^|v) i^j 

and 

p(x.rx.!v) = p(x.!x./v).p(x.|v) = p(x.|v),p(x !v) 

1 J 1 J 3 1 j 



And from these relations# 



P 



x/v 



(x 




p ( X I V ) 




a 

i 



g V ) 
i 



i=l 



i=l 



and finally# 



p ( X ) 
X — 




'' i=l 



6 ( X . 



1 



a . - g . V ) .dv 
1 1 



(A. 12) 
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APPENDIX B : COMPONENTS OF COVARIANCE MATRIX 



For a single measurement of the form z = H,x + Vf the 
observation matrix becomes the vector 




where n is the dimension of the system. 
The Kalman gain is given by 



where P' is used to represent the prediction error covari- 
ance matrixf which is symmetric^ and r is the standard devi- 
ation of the measurement error. 



T T 2 

G = P'H (HP'H + r 1 



(B.l) 



T 

Since H is in this case a vector, the product HP'H 



will be a scalar. 



Let 



c = (HP'H^ + r^r^ 



(B.2) 



The vector P'H can be shown to be 




(B.3) 




and thus the constant c is given by 
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c = 



(B.a) 



n n 2 

1/(E EP'hh + p] 

i=l j=l ij i j 



From (8.1) one now sees that the Kalman gain is a vec- 
tor whose comoonents are the components of P ' multiplied 
by the constant Cr i.e. the components of the gain vector 
are 



n n n 9 

g = ,E p', h. / [ .E ,E P.'.h.h, + r^J 

8 J =1 mj J 1=1 J =1 ij 1 J 



(8,5) 



The estimation error covariance matrix is computed from 
the equation 



P = (I - G.H] .P* = P' + AP 



(B.6) 



where 



AP* = - GHP* = - cP'H HP 



(B.7) 



Let's define 



n 

a. = .E, p' h. 
1 J=1 ij 2 



= o ' h + o.' h + 
il 1 12 2 



. . + p I h 

in n 



(B.8) 



TheOf f rom (B. 3) / 



HP' = (P'H^J ^ = (a, a-, ... a J 

1 z n 



(B.9) 



and from (B.^)/ 



n 



c = I.E, h, a . + r J 
1=1 1 1 



2,-1 



(B.IO) 



From (B.7) the increment aP' to the error covariance 
matrix is t hen 
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individual terms of the estimat 
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APPENDIX C : CHARACTERISTICS OF MEASUREMENT FUNCTIONS 



in Figure C-1 some curves of the form f(jc) = constant 

Si 

are shown, where f(x.) is assumed given. The vector £ and a 
constant value d are also given. The vector is defined 

by the equation 



X 



b _ 





(C-l ) 



or 

x^ = x^ h , Id/ ! h ] 



(C-2) 



where h^ is the gradient of f with respect to x taken at 
point (j<^ , x_^) , and h/|h} is a unit vector along h. 



ill ” ^ o e 



2-2 






Ihj = [h^ + 






lo determine the location of numerical data is 
required. However, it is possible to determine qualitatively 
whether x^ ends on the curve f(£^) t d, or on one side or 
the other of this curve. This is done in the following 
paragraphs for three special functions. 

Case ( i ) - For linear f(x), x*^ always ends on the curve 
fCj(I =flx^) t d, as shown below 
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Figure C - 1 : Problem geometry. 
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Let t(j^) = px^ + qx^ 

I hen = p , h^ = q and jhp = p^ + q^ 

h rom Equation (C-2)» 
x.^ = "*■ qe^J.d/Cp^ + q^ ) 

so 

x^ = x^ t pd/(p^ + q^) 

X 2 = * 2 * qd/(p^ + q^) 

Now # 

f ) = p5^ t qx^ = px^ + p^d/(p^ + q^) + 

+ qx^ t q^d/(p^ + q^ ) = 

= px^ -td = f + d 

Case ( i i ) “ For f(^) defined as belowr x^ always 
on top of f ) + 6t as in the Case (i). 

f(x) = t(x, - p)2 f (x, - q)2]l/2 

— 1 2 

Now / 

h^ = ( x^ - p)/m = ( x^ - q)/m 

where 

m = Kx^ - p)2 

and 

J h;2 = ^2 t h^ 



+ ( x^ - f(x^) 



= 1 



ends 
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= _x^ t I ( ^ ^ *2 ” .d/m 



so 



b _ a 

Xi - x^ 


+ hj^d = x^ ^ ^ 


- p)d/m 




b . a 

X * X 

2 2 


^ h^d = ^ 


- q ) d/m 
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+ ()? t h d - 
2 2 


q,2)V2 = 
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2hid(x^ - p) f 
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Case ( i i 


i) " For f(x) as 


defined below 


and shown 


Figure C-2/ 
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[ d j < TT » 

f(j^) = arctanlCx^ - q)/(x^ - p)] 
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i n 
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Figure C - 2 ; Geometry for Case (iii) . 



f (x^) + d 
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I he case |d| > tt/2 is not of interest since the lines h 
and f(x_J = f ) + d are diverging. For )dj < -n /Z, 

h =“ (x^ -q)/n h = (x^-p)/n 

12 2 1 

where 

n = (xf - p)^ + (xf - q)^ 

• 1 -2 

and 

|h|2 = h^ + = 1/n 

From Equation (C-2)» 

x^ = x^ + d.(”(x^ - q)e + (x^ • p)e ] 

- - 2^-11 -2 

so t 



Now # 



f ( x^ 


) = 


■L 

arctan [ (x° 

2 


- q)/(xb - p)] = 








= arctan 


t(x^ - q) ♦ (x^ - 
2 1 


p)dJ / 








C(x^ - d) - (x^ 
1 2 

4 


- q ) d] 


or 










t(/ 


) = 


arc t an [ (s 


+ d)/(l - sd)J 




where 






• 




s = 


-5 


- q)/(x^ - 


p ) = t an 1 f ( x^^ ) ] 





= X^ - (xj - q)d 
= x^ + (x^ - p)d 
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so 9 



tan[t(x.*^)) = (s + d)/(l - sd) = 

= ltanlf(x.^)J + dJ / tl - d.tan [f (£®) J J = 
= tan[f(xf) + arctan(d)) 
or 

f J = f (jt^ ) t arctan(d) 

Since |d) < ir/2» then arctan(d) < d and ^ will 
reach the line f ( x + d. 



not 
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APPtNUlX D ; ANALYTICAL EVALUATION OF KALMAN FILTERS 



After the design of a filter an analysis phase is 
necessary to verify its behavior in various representative 
si tuat ions. 

A "filtering situation" is considered here to be com- 
pletely specified by: 

--- the true oarameters and initial conditions of the 
system which generate a unique track x(0)/ x(l)f .../ x(t ). 

--- the true parameters of the sensors and their meas- 
urement schedules, 

--- the parameters and initial conditions assumed by 
the filter for the system and for the sensors. 

Ihe approach normally used to determine filter behavior 
is to simulate the desired situation and execute hundreds or 
thousands of Monte Carlo runs and compute sample statistics. 
The most useful results of this process are: 

^™(,kj - the sample mean of the estimated state vector 

at time tj^ obtained from m Monte Carlo runs, 

= ^™(k) - j< ( k ) - the sample mean of the 

estimated error vector at time tj^. 

^™lkj - the sample second central moment of the state 

(or error) estimates about a™ ( k ) , at t, . 

— k 
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(he objective of this section is to study the possibil- 
ity of obtaining values of a^( k ) , _e(k) and F(k)^ which are 
approximately the values of ^™(k)f ^™(k) and ( k ) when m is 
very larger without the use of the time consuming Monte 
Carlo operations. 

1 --- Theoretical Solution 

Unce a "filtering situation"^ as stated above/ is 
defined/ the operation of a Kalman Filter is described by 
the recursive equations below where only the dependence of 
each term on the previous estimates is stressed, 

i(k + l Ik) = f.(£(k) ) tec ! J 

H(ktllk) = (f.(i(k)).P(k),<|)'e(^(k)) t 

t g(£(k) ) ,Q(k) ,g'^(£(k) ) tO.2) 

GU (ktl I k) /P(k + 1 ; k) ) = P(k+1 I k) U (ktl I k) ) . 

lH(i(k>l I k) ) ,P(k+l ; k) .H^(i(k+l I k) ) + R(k)]“^ 

(D.3) 

£(k + l) = ^(k + l!k) + G(£(k+1 ! k) /P(k + p k) ) . (ktin t 

+ v(k + l) - h(x(k + l Ik))] (D.il) 
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K(k+l) = [I - G(^(k+l!k),P(k+l!kJ).HC^(k+llk))].P(k+i;k) 

(D.5) 

in this case these equations describe a deterministic 
process were it not for the random measurement noise v(k+l)/ 
the only random input to the "filter dynamics" since a 
unique track is considered. The initial condition of the 
filter^ £l0) and P(0) are also deterministically given. 

tquations (0.4) and (0.5) can be rewritten as 
i(ktl) = <j) (£(k + l ! k),P(k + l ! k) ) ♦ (£(k + llk),P(k + l!k)).£(ktl) 

P(k + n = f (£(k + l 1 k) rP(ktl ! k) ) 

where (P f f and r are generally nonlinear matrix functions 
of the variables indicated. 

Considering (0.1) and (0.2) one can then write 

i(k+l) = f_* ti^k) ,P(k) ) + g_*(£(k) ,P(k) ) .v_(k+l) (0.6) 

P(k + n = (£(k) ,P(k) ) (0.7) 

^ ^ 

‘where f_» also general ly nonlinear matrix func- 

tions. From (0.6) and (0.7) one can clearly see now thatf 
if £ is a discrete white Gaussian noise procesSf the joint 
process (_x»P> is Markov. 

If one considers a new vector £ formed with the n ele- 
ments of X and the n,(n+l)/2 distinct elements of the (nxn) 
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symmetric matrix Pr y is an tn . ( n + 3 ) /2] -d i mens i ona 1 process. 
Equations (D.6) and (D.7) can thenbe combined into 



where A is a vector function and B a R to R™ matrix func- 
tion/ m = n,(n+3)/2. 

The complete behavior of the filter can be described by 
the probability law of which is obtained from 



and p k + 1 ) J y ( k 1 ) is obtained from p(v(k+l)) and Equation 
(D.b) in the manner shown in Appendix A, The initial value 
^(Oj is deterministic since it contains the initial condi- 
tion of the filter which is given, 

2 --- The Linear Case 

For the special case of linear dynamics and observa- 
tions/ the solution can be obtained in a simpler way. In 
this case the functions (f> / £ and H are not functions of the 
state estimates and so neither is G, 



y(k+l) = A(y(k)) ♦ B ( y ( k ) ) , V ( k+ I ) 



(D.8) 




(D.9) 



Equations (D.l) - (0,5) can be grouped now into 



x(k + l) =<j>(x(k + l!k)/P(k + l/ k)) + 4'(P(k + llk)).v(k + l) 



P(k+l ) = r (P(k+1 ! k) ) 



and finally/ in a recursive way/ 
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(D. 10) 



x_(ktl) = ^*(^(k) ,P(k) ) f g* (P(k) ) .^(k) 

P(ktl ) = h*(P(k) ) (0.11) 

From (O.ll) it is clear that the estimation error 
covariance matrix is now a deterministic process and can be 
precomputedf together with the gains. 

tquation (D.lO) then becomes 

i(k + l) = f_*(i(k)) + a* (k+l) .v(k + l) 

From (D.l) - (0.5) the values of _f and £* can be found 
to be# for this linear case# 

iUtl) = [I f G(k + l).H(k+D) .(()(k).^(k) + 

f G(k + 1 ) .H(k + 1 ) .J(ktl ) + G(k + 1 ) .v(ktl ) 

or 

x_(k+n = S(k + l).x.(k) + F(k+1) + G(k + l).v(k + l) (0.12) 

where 

S(k) = (I + G(k) .H(k)J .<|) (k) 

F(k) = G ( k ) . H ( k ) .j<^( k ) 

the sought after behavior of the filter can be 

described by the moments below# which agree with (121 
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a(k+p = E[^(k+1)] = S(kf 1 ) .E (£(k) J t F(k+1) + 

+G(ktl).E(v(k+l)] (D.13) 

e(k + n = £(kfl) - ^(k + l) = S(k + l).Et£(k)l + 

+ [I - G(k+1 ) ,H(k+l )] .^(k+1 ) + 

+ G(k + l).Etv(k + U) = 

= S(k+l).l(k) + D(k+1 ) .^(k+1 ) - S(k+l).£(k) ♦ 

+ G(k + 1 ) ,E [v(k + l )] (D.14) 

b(k + n = S(k + n .b(k) .s'^Ck + l ) + 

f G(k + l).El^(k + U.vT(k + l)]G^(k + l) (D.15) 

i General Case < 

For the general nonlinear problem one returns to Equa- 
tion (D.9J which cannot normally be solved in closed form. 

Numerical techniques and approximations would be neces- 
sary that may use more computing power anq present worse 
results than the Monte ‘Carlo process that we are trying to 
avo i d . 

One way to simplify the problem is to assume that y(k) 
has an approximate Gaussian distribution/ even with all the 
non 1 i nea r i t i es shown. Linearizing Equation (0.8) around the 
mean value of ^(k) would give 
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X(k + 1) = A(y(k)) + A(v(lc) ) . [y (k) - y(k)) + 

+B(^(k)).^(k+l) (0.16) 

Taking ^(0) as the deterministically known initial 
value ^(0)f and applying the expectation and variance opera- 
tors to Equation (D,l6)f one gets the recursive equations 
for the moments of ^(k): 

^(k+n = A(7(k)) ♦ B(7(W) ) .E (v(k+l)l (0.17) 

Varly(kfl)J = [ ^ A ("^ ( k ) ) ) . Var [y ( k ) J . ( -i- A(y(k))l'^ + 

+ B(y(k)).Var lv(k + l)J .B^(y(k)) (0.18) 

Ihe moments we are interested in, F, ^ and ^ are 
directly obtained from subvectors and submatrices of the 
above moments. 

The primary difficulty with this approach, however, is 
in obtaining the functions A and B. This can be seen for the 
very simple case below, as an example. 

Suppose a scalar system with a single observation is 
described by the equations 

x(k+l) = x^(k) + w(k) 

z(k) = x(k) ♦ v(k) 



18 a 



I 

I 

I 

I 

I 

« 






I 



tquations (0.1) - (0.5) give 

i(k+l ! k) = i 2(k) ; <j)(k) = 2x(k) 

p(k + lJk) = a;^(k).p(k) + q^(k) 

(k) .p(k) + q^(k) 

0(k + l) = 5 r 

" ■ ' ai^(k).p(k) + q^(k) + r^(k) 



and/ after the appropriate steps 



X (kf 1 ) 



p(k + l ) 



X ( k ) + 



yx^(k).p(k) + q^(k) 

4x^(k) .p(k) + q^(k) + r ^(k) 



. (x(k) - x(k)] 



(4i^ (k) .p(k) + q^(k)l . r^ (k) 

ai^(k).p(k) t q^(k) f- r^(k) 



ai^(k) .p(k) + q^U) 

4x ^(k) .p(k) + q^(k) + r^ (k) 



.v(k + l ) 



fhe mean value given by Equation (0.17) can be obtained 
in a simpler way/ however. It is the result obtained by 
running the filter in the simulated environment but making 
the measurement noise assume its mean value E(v(k))/ nor- 
mally zero. 

If one uses more terms in the expansion of Equation 
(0.6)/ better results will % be obtained at the expense of 
increased computing effort and time. 
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For this general case it is suggested that these equa- 



tions be applied to a simple scalar or two 
linear problem and that the computing power 
the sought after moments be compared with 
run a Monte Carlo process yielding the same 
racy in the results. 



dimensional non- 
requ i red to find 
that requi red to 
level of accu- 
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APPENDIX E ; COMPUTER PROGRAM 



Ihe basic structure of the computer program written for 
the simulation of the tracking problem^ and the evaluation 
of models and filtering algorithms^ is shown in Figure E~l. 
The solid lines represent the normal flow of control within 
the program; the dashed lines show the extra pathes that 
become available whenever an interruption is requested by 
the user. 

After a brief introduction to the program a table 
called MENU is presented to the user/ as shown in Figure E- 
2. if the user chooses actions number 1 or 5 the result is 
clear. 

Choice number 3 provides access to all the problem 
variables. Since the number of variables that characterize 
each simulation is very large/ a simple way to deal with 
these variables had to be devised. As it is implemented/ 
the program always "remember" the values given to the vari- 
ables at the last time the program was used/ so that only 
the variables to have their values modified have to be 
addressed. The change of value of a variable is simply made 
by selecting the page where it appears/ typing the letter 
associated with the-variable and the new desired value. The 
program immediately responds by presenting back the entered 
value. A more detailed diagram of the actions resulting 
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F i au re E - 1 



Basic srpucture of oroaran^. 
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MFfJll 



Press the nuTber corresoonding to the action 
des i red and c / r . 

t - PRFSENT tutorial I^!FnR^>iAT ION 

? - ^>^O^IFY program flags 

3 - formulate or modtfy the problem 

U - START THE SOIUTIOM OP THE PROBLEM 
5 - FNO the program A.mO «^XTT 



Fiaure E-B J IFe MENU 



189 



from choice number 3 is presented in Figure E-3. A typical 
page of variables is shown in Figure E-^. 

Choice number ^ simulates the tracking problem accord- 
ing to the values defined in the previous step. A block 
diagram of the actions involved in the simulation is 
presented in Figures E-5/ E-6 and E-7, 

Initially all the important problem variables are 
printed with their current values in a form as shown in Fig- 
ures E-6/ E-9 and E-10, During the first run the scheduling 
of all the events involved in the simulation is also printed 
for future reference/ as shown in Figure E-ll. 

With choice number 2 from the MENU a new table is 
presented as shown in Figure E-12. This table in also 
presented to the user whenever he request s an interruption/ 
at any time or point within the program. 

The extra choices that now become available are almost 
self-explanatory but it should be added that with choice 
number 7 the simulation is ended and the partial statistics 
are computed and written in the appropriate output files; 
choice number 6 allows the modification of any variable in 
the middle of the simulation/ in the same way as explained 
before/ choice number 0 transfer control back to MENU or to 
the point of i nterrupt i on . 

A listing of the program/ in C language/ is available 
at the Electrical Engineering Department/ Naval Postgraduate 
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to a to 

CPC ’ MEMU 



FiaureE”3 t Chanqino the value of problem variables 
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3 






To change the value of any variable/ type the indicated 
letter/ at least one scace/ the desired value and c/r. 

To see the other variables of the Drobleip/ tvpe 1 c/r. 

To qo back to menu, tyoe 0 c/r. 

( 1 meter = l.hqy yards; 1 meter/sec = knots ) 

TARGFT IMITIAL PARAMfTFRS 



initial x - no?^ition 


: 0.0 


km 


a 


y • Dosition 


: 0.0 


km 


b 


speed 


: 7.0 


m/sec 


c 


head i no 


: 30S.0 


deq 


d 


frequency 


: SOO.ft 


hertz 


e 



standard deviation of 


f o r c i no 


functions 




speed/sec 


: 0.0 


m/sec/sec 


f 


headi no/sec 


: O.n 


deq/sec 


q 


f reouencv/sec 


T 0.0 


hertz/sec 


h 



Figure l Tvoical cage of variables. 
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from MFNlj 







Figure E -5 



Simulation T 
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Fioure E-b ; Simulation II 
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Fiaure E-7 : Simulation Til 
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TARGfT IMITIAL PAHAKFTFRS 



initial y( - oos i f i on 


0.0 


K m 


a 


y - no?tition 


0.0 


kn 


b 


speed 


7.0 


m/sec 


c 


head i nq 


■^os.o 


deq 


d 


f requenc V 


=^00. fl 


hertz 


e 


standard deviation of 1 
soeed/ sec 


' 0 r c i no 
0.0 


f unc t ions 
m/ sec /sec 


f 


headi no/sec 


0.0 


deq/sec 


o 


f requenc v/sec 


0.0 


hertz/sec 


h 



GENFRAL niRECTIOMS 



nu'TDe r 


of 


n I F AR buoy s. ! 


r 


buoys 


numbe r 


0 f 


LOFAP bucys : 


: 0 


buovs 


number 


of 


maneuve r s 


? 


man 


number 


0 f 


runs 


?o 


runs 


durat i on 


of runs 


2000 


sec 


numbe r 


o f 


Icey points 


20 


points 


seed 
a ve r aqe 


sound ve 1 oc i t y 


1 500.0 


m/sec 


range of 


the buovs 


8.0 


km 



filter initial PARAMFTFRS 



initial x • nos i t i on 


-0.5 


km 


a 


y - nosition 


n .0 


k m 


b 


«^oeed 


6.0 


m / s 


c 


^ e a a i no 


^ ! 0 . n 


d ;- 


d 


f renuencv 


5 0 0.5 


h e r t Z 


e 


standard deviation of fore i no 
spped/sec ! O.OOS 


f unc t i on s 
m/sec /sec 


f 


heaainn/ser ! 


: 0.050 


oeq/sec 


n 


f requenc y/sec : 


: 0.0 10 


hertz/sec 


h 



Fiaure fc-W : LisTinq of oroblen» variables I 
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(X) Jj ox 



initial COVAPIAi\CE matrix 





X DOS 


V nos 


soeed 


heaoa 


f r en 


X pos 


1.0 a 


0.3 b 


0.0 c 


0.0 d 


0.0 


y DOS 


* 


1 .0 q 


0.0 h 


0.0 i 


0.0 


speed 


★ 


* 


3.0 m 


0.0 n 


0.0 


heada 


* 


* 


* 


100. 0 s 


0.0 


f req 


* 


★ 


* 


* 


1.0 


(km/ 


m/ sec / 


rieorees and 


he r t z 


c rossmu 1 t i o 1 


i ed ) 



measurements schedule i 



buoy 


t y oe 


start 


period 


mean 


St dev 




1 


B 


100 


1 00 


0.0 


5.0 


abc de f 


1 


F 


mo 


100 


0.0 


0 .oa 


nh 1 j k ) 


? 


e 


500 


1 00 


0.0 


5.0 


mnoDo r 


? 


F 


500 


1 0 0 


0.0 


o.oa 


s t u V w X 


(start 


and 


period 


i n seconds ) 







(mean and std dev in degrees# hert? or seconds) 
(period must be ^ero i^ measurement is not used) 



MEASUREMENTS SCHEDULE II 



buoy 


type 


start 


period 


mean 


St dev 




1? 


T 


500 


100 


0.0 


0.01 


abcde t 


IS 


T 


500 


1 00 


0.0 


0.01 


o h i j k 1 


33 


T 


500 


1 00 


0.0 


0.01 


m n o p o r 


3 


B 


500 


1 0 0 


0.0 


5.0 


S t U V W X 


(start 


and 


pe r i od 


in seconds) 







(mean and std dev in degrees# hertz or seconds) 
(period must oe zero if measurement is not used) 



Eioure E-9 i Listina of nrnblem variables II ' 
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TARGET MAfMFUVEPG 



1st maneuver : time = 600 sec a 

vriot = 0,0 m/sec/min b 

hdot = -R.O deq/min c 

fdot = 0,0 hertz/min d 

TARGET MANEUVERS 

?nd maneuver : time = 1?00 sec a 

vdot = 0,0 m/sec/min b 

hdot = 0,0 deq/min c 

fdot = 0.0 hertz/min d 

RUOYS parameters 

Eirst buoy : tvpe = DIEAR a 

x-position= 0.0 Vfr b 

y-position= -3.0 c 

If DIFAR: 

bearing error mean= 0,0 deq d 

std dev = S.O aeg e 

BUOYS parameters 

(?nd buoy : tvoe = DIPAP a 

x-position= '4.0 Wm b 

y-position= -5.0 km c 

If DTFAR: 

bearinqerrormean= 0,0 oeg d 

sto oev = S,0 cjoq ^ 

Eioure E-10 : List i no of oroblem variables III 
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SEQUENCE OF EVENTS FRUV 601 TO 900 SECONDS 



700 


sec 


— 


bear i no measurement by 


buoy number 1 




700 


sec 


- 


frequencv measurement 


by buoy number 


1 


700 


sec 


- 


bear i no measurement by 


buov number 2 




700 


sec 


- 


frequencv measurement 


by buoy number 


2 


700 


sec 


- 


time delay measurement 


bv buoys 1 and 


? 


700 


sec 


- 


^^onte Carlo roint number 7 




«00 


sec 


- 


bearino measurement by 


buov number 1 




flOO 


sec 


- 


frequency measurement 


by buoy number 


1 


FOO 


sec 




bearino measurement by 


buoy number 2 




FOO 


sec 


- 


frequency measurement 


by buov number 


2 


BOO 


sec 


- 


time delay measurement 


by buoys 1 and 


2 


BOO 


sec 


- 


^onte Carlo coint numb 


e r 7 




900 


sec 


- 


bearino measurement by 


buov numbe r 1 




900 


sec 


- 


frequency measurement 


by buov number 


1 


900 


sec 


- 


bearino measurement by 


buoy number 2 




OQO 


sec 


- 


frecuency measurement 


by buov num.ber 


? 


900 


sec 


- 


time delay measurement 


by buoys 1 and 


? 


900 


sec 




Monte Carlo COint number 7 





Fiaure E-11 : Table of ever'ts. 
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PROCPAN- CHARACTFPISTICS 



Press the nurnber corresnonHinq to the option or 
action desired and c/r. 

To continue or return oress 0 c/r. 



1 • update tarcet everv second 



? - uodate taroet only before measurements 



3 •" start orintinq on line or inter 



^ - stop print inn 



5 - out results of next run on outout files 



6 - modify parameters of the problem 



7 - terminate the problem 



Finure E.-l? t Interruption table* 
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School t 



to any interested reader. 
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